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Abstract 


The study of ultraluminous X-ray sources (ULXs) has changed dramatically over the last decade. In this review we 
first describe the most important observations of ULXs in various wavebands, and across multiple scales in space and 
time. We discuss recent progress and current unanswered questions. We consider the range of current theories of ULX 
properties in the light of this observational progress. Applying these models to neutron-star ULXs offers particularly 
stringent tests, as this is the unique case where the mass of the accretor is effectively fixed. 


1. Introduction 


Ultraluminous X-ray sources (ULXs) are usually defined by fluxes which, if assumed isotropic, give luminosities 
Lx > 10? ergs |. (1) 


The implication is that ULXs appear to be above the Eddington critical luminosity 


M 
Leda = 1.3 x 1078 (5) ergs |, (2) 


© 


for accretor masses M < 10 Mo (see Eq. [20|for the exact definition of the luminosity Lgaa). Early systematic X-ray 
surveys found significant numbers of these objects: reported 16 sources with Ly > 10% ergs~!. By 
now, more than about 1800 ULXs 2022) are known, including several with Ly > 10*! — 10” erg s™!. 
ULXs are (almost) all in external galaxies, but not located in their nuclei, making them distinct from accreting super- 
massive black holes 

From the beginning it was clear that the brightest ULXs required a significant shift in the standard paradigm for 
the then known accretion—powered sources, and this provoked a wide variety of speculations. An obvious possibility 
was that the accreting mass M was unusually large. (Colbert and Mushotzky] {1995} used a sample of ~ 20 ULXs with 
apparent X-ray luminosities in the range 10*° — 10° ergs“! to argue that the typical accretors were black holes with 
masses in the range 107 Mo < M < 10*Mo. The suggested objects soon acquired the name ‘intermediate—mass black 


holes’, or IMBH (Taniguchi et al.}|2000). Models for IMBH||formation include rapid merging of stellar—-mass black 
holes in globular clusters (Miller and Hamilton} |2002), or of stars in dense clusters (Portegies Zwart and McMillan 


'Some dwarf galaxies have active nuclei powered by accretion on to black holes of masses ~ 10° Mo, and so with similar luminosities to some 
ULXs. These black holes appear to obey the same M — o relation between black hole mass and host spheroid velocity dispersion as universally 
holds for higher mass nuclear black holes and their host galaxies (see[Baldassare et al [2020][King and Nealon]2019]for a discussion). Unlike these 
objects, ULXs are not unique systems at the dynamical centres of their hosts. Accordingly, in this review we confine use of the term ‘IMBH’ to 
objects not in the nuclei of dwarf galaxies (see Section[2.7}. 
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2002). A problem for all these models was to ensure that the IMBH would find a companion star and then form a 
binary able to sustain mass transfer at the required rate (up to ~ 1075 — 1074 Mo yr7!). 

A quite different idea was that ULXs were fairly standard X-ray binaries with M < 10M in unusual states. 
suggested that the accretion discs in standard black hole X-ray binaries might be unstable to the 
photon-bubble instability, and simply radiate at super-Eddington luminosities. This is limited to luminosities Lsph < 
3 x 10% ergs!, and requires that mass loss from the disc surface is small. suggested that 
the high luminosities of ULXs could result from relativistic beaming in jets from microquasars oriented towards the 
observer. This is again limited, to Lsph S 10% ergs"!. Similarly, it was known that magnetic fields with the strength 
~ 10! — 10! G found in accreting neutron stars reduced the electron scattering opacity for X-rays propagating 


perpendicular to the field 19717 but this works only for Lenn S 3 X 10° erg s~!, and even then requires 
a magnetar-strength field B ~ 10'°G (Mushtukov et al. [2015b), although these do not seem to appear in binary 
systems (see Section|3.9). The limits for all three of these stellar-mass models are incompatible with the luminosities 
Lx > 104! — 10%? erg s~! observed in some ULXs. suggested that radiatively inefficient accretion 
on to black holes could explain ULXs as the radiative output was tightly beamed. Here much of the accreting gas is 
swallowed by the hole without radiating. 


(2001) noted that the non-magnetic neutron star in the low—mass X-ray binary Cyg X—2 (not a ULX) 
was known from evolutionary arguments (King and Ritter} |1999} |Podsiadlowski and Rappaport! |2000} see Section 
0 


3.11] below) to have survived a previous phase where it was subject to mass transfer at rates ~ 107 — 10° times its 
Eddington accretion rate of ~ 1078 Mo yr™!. The neutron star had evidently ejected the vast bulk (~ 2 — 3 Mo) of the 
transferred mass, retaining no more than a few x 0.1 Mo. Ejection like this is expected in the simple picture of 
(1973) of mass expulsion by radiation pressure from thin accretion discs fed at super—Eddington rates”} 
(see Section below), so a plausible deduction was that the intense accretion disc wind might block the emitted 
X-rays from many lines of sight, perhaps leaving only narrow escape routes near the disc’s rotational poles. If these 
had total solid angle 4b, with b « 1, an observer viewing the source along these directions but assuming that the 
emission was isotropic would ascribe a luminosity 


1 
Lei = ma >L (3) 


to the source, where L is its true luminosity] 
Early discussion of this model was hampered by the lack of a clear relation between the beaming factor b and 
other physical variables, which left a spurious degree of freedom. This gap was later closed (see Section |3.6), but 


even before this the model made interesting and testable predictions. Simple estimates showed that the characteristic 
velocity of the accretion disc wind was ~ 0.1c (King and Pounds] 2003| compare Section (3.14). If b < 1072, the 
X-ray fluxes of the early observed ULXs gave L = bL spn < 10°’ erg s~’, below the Eddington value for a 10 Mo black 
hole. King et al (2001) also noted that for b < 107°, the deduced L would be compatible with a neutron star accretor, 
as expected if Cyg X-2 is a survivor of a ULX phase, and with a white dwarf accretor in cases where the observed 
X-ray flux was soft (photon energies < 0.1 keV). The brightest ULXs would correspond to the same value b = 107° 
of the beaming factor, but with a 10 Mo black hole accretor. 

This idea gave sensible answers for the likely populations and lifetimes of ULXs (King et al.|[2001). It suggested 
that ULXs were an extremely common stage of high—mass X-ray binary evolution characterised by very large mass 
transfer rates. 

From the analysis of Cyg X-2’s evolution, the best candidate for this stage was the preceding phase of super- 
Eddington mass transfer. This inevitably follows the standard wind—fed high—mass X-ray binary stage once the donor 
star fills its Roche lobe, as this is always more massive than the accretor. Then mass transfer shrinks the binary to 
conserve orbital angular momentum, maintaining very high mass rates limited only by the companion’s ability to 
expand rapidly, usually on its thermal timescale (see Section [3.1 1p. Begelman et al. (2006) noted that the extreme 


>The possibility of this effect is probably the reason why the occasional Ly ~ 10°? ergs~! outbursts of the Be-star — magnetic neutron star 


system A0538-66 did not attract more attention. 

3Note that shes Sar Sma 73} cannot change the accretion luminosity for a neutron star accretor. 

4 Although this is almost always called ‘beaming’, as in this review, it is important to note that it is simply collimation, rather than more complex 
processes such as relativistic beaming. 


binary SS433 was probably a system of this type, but not viewed along the accretion disc axis, so was probably a ULX 
‘viewed from the side’. 

The current value of the mass transfer rate My in any system like this is (as usual) hard to determine from direct 
observations. In particular measuring the rate P of change of the orbital period P gives at best a wide upper limit 
My < |È/P|M; to the mass transfer rate (where M3 is the mass of the donor star). This is clearly seen for high—mass 
X-ray binaries: measure P/P ~ —10~° yr"! for LMC X-4, Cen X-3, 4U1538-522 and SMC 
X-1, which would give maximum mass transfer rates ~ 1076 Mo yr! and make all of these systems ULXs, but their 
observed luminosities are all sub—Eddington (My < 1078 Mo yr7!). 

We now know that a small group of high-mass X-ray binaries with Be—star companions can occasionally produce 
very high mass transfer rates for short intervals because of dynamical effects on the Be-star disc — see Section [3.15] 
These systems are unique in repeatedly making transitions from a ‘normal’ X-ray binary to a ULX and back. This 
property distinguishes very sharply between models for ULXs. 

The short lifetimes of both types of high-mass X-ray binary (pre-ULX, or Be-star) showed that the disc-wind 
beaming model required an association of ULXs with star formation. This was already suspected, in particular fol- 
lowing the discovery of 7 ULXs in Chandra observations of the Antennae, where the merger of two galaxies has 
provoked vigorous star formation (Fabiano et al.| 2001). Chandra observations of the Cartwheel galaxy 
soon strongly confirmed it. They revealed more than 20 ULXs with Lsph > 3 X 10% erg s~}, four of them in the 
range 1—5.5x 10*! erg s™!, in a spreading ring of star formation triggered some 3 x 108 yr ago when a smaller intruder 
galaxy crossed the plane of the Cartwheel disc close to its nucleus. None of these features were easily compatible 
with the IMBH hypothesis for ULXs (King||2004). The proof that at least some ULXs had stellar masses, and indeed 


Cartwheel (X-ray contours on HST/WFPC2) 
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Figure 1 Broadband X-ray contours overlaid on the HST/WFPC2 optical image. The lowest contours are 0.0345, 
0.0431 counts pixel”! (pixel=~ 0”.5) and then increase successively by a factor of 2.(Gao et al. 2003). 


neutron star accretors, as suggested by (2001) (see the discussion after Eq. (3), finally arrived when|Bachetti| 
(2014) found a coherent periodicity P = 1.37s in the ULX M82 X-2, with Lyn = 1.8 x 10% ergs~!. The only 


plausible interpretation, as the spin of a magnetic neutron star, required M ~ 1 Mo, and showed that this object had 


3 


an apparent X-ray luminosity Lspn more than 100 times its Eddington luminosity. Note that this observation did 
not identify the cause of the neutron star’s super—Eddington luminosity — this is still disputed, essentially between 
geometrical beaming as suggested by (2001) and the effects of very strong magnetic fields — see the later 
part of the present review, and the reviews by and/Mushtukov and Tsygankov|(2022). Models of 
this type suggest that a significant, if not dominant, fraction of ULXs are accreting NSs (see abstract in/Mushtukov| 
fet al.[2015a). It is interesting that only 6 out of the ~ 1800 non—Be ULXs have so far been observed to show pulsing, 
unlike what we might expect if this dominance holds. (We thank a referee for provoking this insight). A further 
problem here is that these very strong fields do not seem to be present in binaries which do not have high accretion 
rates. This picture then requires us to believe that strong—field neutron stars in binaries are unobservable unless the 
binary companion supplies them with mass, at a super—Eddington rate. 

The observation by [Bachetti et al.|(2014) removed the main argument motivating IMBH models, and discussion 
of ULXs is now tightly focussed. The finding of ULXs with Ly > 104! — 10” ergs™!, at least one of which contains 
a neutron-star accretor is a stringent test of theoretical models. It is compatible with disc—wind beaming, and this is 
also the only model naturally allowing standard Be — X-ray binaries to become ULXs during bright outbursts. 

The main question in discussions of ULXs is whether they are all generically similar, or divide into subgroups 
with different physical models for each. We discuss this question in detail at the end of this review (Sect. ø. but note 
here that in the absence of observations dividing ULXs into clearly distinct groups, Occam’s razor |leaves disc—wind 
beaming as the only model open to test against all observations of ULXs. For reviews taking differing points of view, 


see|Fabrika et al.|(2021), and|Mushtukov and Tsygankov| (2022), the latter for pulsed ULXs containing very strongly 


magnetic neutron star accretors. 


2. Observations 


Following the launch of NASA’s Chandra (Weisskopf et al.|2000) and ESA’s XMM-Newton (Jansen et al./2001) 


satellites in 1999, the field of studying X-ray binaries in external galaxies was revolutionised. Both satellites carry 
X-ray CCD instruments sensitive in the 0.3-10 keV range, and have spectral, timing and imaging capabilities (with 
on-axis resolutions of < 5”), the latter being vital for successfully disentangling an X-ray bright object from the centre 
of its host galaxy. 


2.1. Catalogues & X-ray Surveys 


After more than two decades of pointed observations, the number of individual sources now totals over 300,000 


from Chandra (CSC 2.0: 2010) and over 500,000 from XMM-Newton (4XMM-DR39: [Webb et al./2020) 


respectively. These extensive catalogues (and their previous iterations) have led to the formation of valuable ULX 


catalogues, most recently by|Kovlakas et al.|(2020) utilising Chandra, (2019) utilising XMM-Newton 
and (Walton et al.|/2022) using the combination of Chandra, XMM-Newton and Swift (with previous catalogues by 


Walton et al. 201T}{Wang et al.[2016). 

In these catalogues, ULXs ( numbering ~ 1800 in the most recent: are located by cross- 
matching against local galaxy catalogues (e.g. HECATE: |Kovlakas et al./2021), with checks that sources lie within 
the host’s optical D25 isophote but are not coincident with the central region, and estimates applied for foreground and 
background (AGN) contamination. These catalogues enable quantification and searches for variability (e.g. propellor 
states), targeted follow-up observations, and the study of correlations with host type, including star formation rate and 


metallicity, thereby testing the predictions of binary population synthesis models (e.g. 2019 
Kuranov et al.[2020). 


Both of the most recent ULX catalogues find that spirals are the most common host galaxies — consistent with 
ULXs being associated with regions of high star formation. An excess of ULXs is found in lower metallicity galaxies 


(see also|Soria et al./2005}|Mapelli et al./2010}|/Prestwich et al.)2013}|Brorby et al./2014} 2016), the 
numbers scaling with star formation rate as (Kovlakas et al.!2020): 


5‘Do not multiply hypotheses beyond necessity’, or more simply, ‘don’t invent two theories for the same thing.’ 
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Figure 2 From {Mineo et al.] (2012) showing the point source X-ray luminosity functions (XLFs) (normalised by 
star formation rate) created from Chandra observations of M101 (blue), NGC 3079 (red), M51 (cyan), Antennae 


(magenta), NGC 3310 (green). 


Nurx = 0.45120% x SFR m, /year + 3.3238 x Ta (4) 
It is important to note that, regardless of spatial coverage, the long timescale (months or longer) variability of 
ULXs (as discussed below) demands repeat observing to obtain a reliable census of the population. The all-sky 
survey by eROSITA (eRASS -{Cappelluti et al./2011} will provide precisely these repeat scans, eight being taken over 
the space of four years (2020-2024) with a far greater sensitivity than the ROSAT all-sky survey (Voges et al.]1999) 
and coverage from 0.3-10 keV (matching the nominal range of Chandra and XMM-Newton). 
The production of such extensive catalogues provides unique statistical insights, notably the X-ray luminosity 


function (XLF, such as that shown in Figure 2) giving the combined numbers of low— and high—mass X-ray binary 


populations and how these depend on metallicty (Lehmer et al.|/2021). In turn these XLFs allow estimates for the 


impact of accreting binary populations on reionisation and galaxy evolution (Fragos et al./2013}|Kovlakas et al.|2022). 
XLFs can also be directly compared with population synthesis calculations which include the evolution of systems 


through their various phases (including mass transfer, common envelope, SNe and gravitational inspiral). Recent 
attempts include those by where the effect of geometrical beaming (following equation 
55) is included, 2021), where the observational impact of precession is considered (see below), and 
(2020) who obtain the XLF for magnetised neutron stars in the absence of beaming and precession of 
the wind-cone (see section 2.3.3). 


2.2. Spectroscopy 

Chandra and XMM-Newton have made it possible to study ULXs with (relatively) high signal-to-noise and energy 
resolution, and have provided a ‘canonical’ picture of ULX X-ray spectra. The improved data quality from these 
satellites allowed single component spectral models to be excluded, instead indicating the presence of a soft excess 
and signs of a break to a steeper spectral index at the highest energies, where instrumental effective area typically 


diminishes (e.g. 2003} Miller et al./2004b 2005} |Stobbart et al.|2006). 
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Figure 3 From|Pinto et al.|(2017), illustrating the variety of observed spectra of ULXs. Each dataset has been unfolded 
through the instrumental response only (XMM-Newton and/or NuSTAR) and show two spectral states for each ULX 


(one hard and one soft). The fluxes on the y-axis are arbitrary. 


If the soft excess is modelled as a thin disc (set to truncate at the ISCO), the characteristic temperatures are suf- 
ficiently low (0.1-0.2 keV) that IMBHs accreting at significantly sub-Eddington rates were proposed 
[2004alc). But these accretion rates imply an upscattered power-law component extending unbroken to very high ener- 
gies (around the electron’s thermal energy, see e.g. Done et al.|2007). The presence of a break below 10 keV conflicts 
with this (as too do theoretical arguments - see SoD cain the alternative hypothesis of super—critical 
accretion (Shakura and Sunyaev|1973). The first large-scale spectral analysis based on a picture of super-Eddington 
mass supply was performed by (2009). Here the harder emission was created by inverse Compton 
scattering of seed photons (presumed in this case to originate from the soft excess). Since then the field has benefited 
from data from newer X-ray satellites, notably NuSTAR, whose ability to image in the 3-70 keV bandpass (Harrison 
et al.[2013) confirmed the < 10 keV break with high signal-to-noise ratio (e.g. 
Taal aoned for the most complete picture of the X-ray SEDs of ULXs (Figure 3). This wide band-pass has 
led to important advances in the field, including the discovery of the first ULX pulsar (Bachetti et al.]2014). 

Current modelling of the X-ray spectra of ULXs remains to some degree phenomenological (although more phys- 
ically motivated models have been developed: e.g. |Vinokurov et al./2013). Models used to describe the soft compo- 
nent variously include an advection dominated disc (typically piskpBs — with a radial temperature profile T œ r~? with 


p < 0.75, e.g. |Mineshige et al.]1994), or a standard disk blackbody (typically DISKBB — 1984). Each of 
these has been used to describe the wind, inclination and advection-modified emission from the super-critical disc or 


from an accretion curtain (depending on field strength and accretion rate —[Mushtukov et al.[2017). The harder X-ray 
emission (typically > 1 keV) is variously modelled with Compton components (with empirical assumptions), broken 
power-laws and — more recently — the addition of emission from a neutron star accretion column (Walton et al.[2018c). 
Such two- or three-component spectral fits are motivated by spectral timing methods (e.g. the covariance spectrum - 
and pulse-resolved modelling Walton et al.[2018b). Although multiple 
components are often used, simple modelling has encouraged the wide adoption of three simple, empirical categories: 
soft ultraluminous, hard ultraluminous, and broadened disclike (Sutton et al.[2013|- Figure 4). It is important to note 
that these definitions are merely descriptive, they do not imply precise statements about the relative contributions of 
the underlying physical components but provide a useful shorthand for the shape of the spectrum and its evolution 
with time. 

The simple phenomenological modelling described above can be connected to the theoretical picture of super— 
Eddington mass supply. The characteristic temperatures obtained from modelling the spectra can approximate the 
colour-corrected temperatures Toph, Tcsp and Tomax, Of the quasi-blackbody emission from the disc at characteristic 
radii (Poutanen et al./2007): the outer photosphere of the wind (Rph), spherization radius (Ry, — noting that the 


assumed position of this radius in 2007| slightly differs from Eq. 14), and innermost edge of the disc 
(Rmax): 


copy" 
Toph = 0.8 fool (2) mm? keV (5) 
Ew 
Tom = US fon "Im"? keV (6) 
Tomax = lef’ keV (7) 


Associated values for the accretion rate then follow under the assumption of a fixed fraction of accretion energy 
deposited in the wind (ew), colour temperature correction (f¢o1) as well as dynamical values associated with the wind 
(i.e. the wind velocity relative to the Keplerian velocity at rsp), 8, and the cotangent of the wind opening angle, € — 
Poutanen et al.|2007). As an example, typical values of the soft X-ray component’s temperature lie between 0.1 and 
0.4 keV (Middleton et al./2015a); for values of foo. ~ 2 this indicates Eddington-scaled mass supply rates in the range 
~ 20-300 for a 10 Mo black hole and 50-800 for a 1.4 Mo neutron star. 

One must be cautious when modelling ULX spectra, as descriptions are often degenerate (see Section 2.3.4), and 
the presence of a neutron star can have a major impact on the interpretation of characteristic temperatures (see Section 
2.2 for more detail). As an example, strong dipole fields could truncate the thin disc (i.e. before the local Eddington 
limit is reached) which comes with the natural corollary that such a disc emits below the Eddington limit for a neutron 


star mass (see|Middleton et al.}2019a). Additional uncertainty in the spectral fitting (regardless of the nature of the 
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Figure 4 From|Kaaret et al.|(2017), showing the variety of phenomenological states of ULXs (SSUL: supersoft ultra- 
luminous; BD: broadened disc; H/SUL hard and soft ultraluminous) as defined by 2013) as compared to 
the standard canonical soft (thermal) and hard state spectra of XRBs. 


compact object) arises as Compton down-scattering within the wind itself can distort the temperature of the photons 
produced at various radii. It remains likely that the soft component is inclination-dependent, so that a more edge-on 
system appears increasingly softer until the temperature from the face of the wind dominates (and from which we 
expect a luminosity ~ Lgaa to be emitted: [Poutanen et al.[2007| [Zhou et al.[2019} [Qiu and Feng]2021). Together with 
changes in accretion rate, the inclination dependence provides a natural explanation for soft and super-soft (e.g. 
et al.[2016) ULXs which fit well into the simplified unified model of ULX spectra 


et al./2015a}|Soria and Kong}2016}|Urquhart and Soria|2016a 2017). 
A notable absence from ULX spectra (with the exception of Swift JO243.6+6124: |Bykov et al.|2022) see also the 


claim of a feature in NGC 7456 ULX-1: is that of Fe Kg emission, with (2012) 
placing strong constraints on its presence in Ho IX X-1. This absence is remarkable, given the strong fluorescent 
yield of Fe and its typically high abundance. Importantly this may point towards an extremely highly ionized gas 
environment, consistent with an irradiated conical outflow (see also[Medvedev and Fabrika|2010}|Pinto et al.[2021). It 
is quite likely that weak Fe K in reflection would be produced far out in the wind cone where the material is not fully 
ionized, in the disc beyond Rspn provided that this radius is small enough that the outer disc is irradiated by X-rays 
with energies > 7 keV, and beyond the magnetospheric radius (so long as this is > Rsph). 

More elaborate models of emission from slim discs already exist (e.g. [Straub et al.[2013), and models for fitting to 
data may now be obtained directly from 3D radiative MHD simulations for super—Eddington mass supply on to black 
holes (e.g. and magnetised neutron stars (Kawashima and Ohsugal2020). Still more elaborate 
spectral modelling may include self—consistent photoionization. Looking forwards, it is hoped that such models may 
bridge the gap between theory and observation. 


2.2.1. High energy-resolution X-ray spectroscopy 

Phenomenological models can reasonably describe ULX spectra as observed by current instrumentation, although 
clear residuals to these models are found between 0.5 and 2 keV in multiple ULXs ( 
et al./2015c), and seen in both XMM-Newton and Chandra spectra (ruling them out as instrumental artefacts - 
et al.|2006; [Pintore et al. 2015). Initially thought to originate in diffuse star-forming plasma (Strohmayer et al.|2007| 
[Miller et al.[2013), the luminosities are too high for such an interpretation, and an 
association with mass-loaded outflows was instead made (Middleton et al.[2014). Both via high angular resolution 
studies with Chandra (Sutton et al./2015), and by studying the correlated variability of the features with spectral 
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hardness (Middleton et al./2015c 2020b), it was determined that they must be associated with the ULX 


and not with the surrounding nebula (e.g. |Sathyaprakash et al.!2019a). 

Instruments with higher energy resolution (e.g. XMM-Newton’s Reflection Grating Spectrometer, RGS) have 
shown that the CCD-quality features are indeed composed of blueshifted absorption lines and rest-frame emission 
lines (Pinto et al.[2016} [2017|- MEE STOO) Ta N apparently blueshifted emission lines in NGC 5204 X- 
1), confirming the suggestion by (2014). In NGC 1313 X-1, counterpart Fe absorption lines have also 
been found in NuSTAR CCD-quality data with a blueshift consistent with that inferred from the lower energy resolved 
lines (Walton et al.[2016b). The limit implied on the outflow velocity is > 0.2 c, considerably faster than observations 
of thermal (or magneto-centrifugal) winds from sub—Eddington XRBs (e.g. |Ponti et al.//2012), but consistent with 
radiatively driven winds from a super-critical inflow (e.g. 
(2021). 

Ultrafast outflows have now been located in several ULXs via individual searches of RGS spectra (NGC 1313 X-1 
and NGC 5408 X-1: NGC 55 ULX: NGC 247 ULX-1: and 


notably in the transient NS ULXs, NGC 300 ULX-1 (Kosec et al.!2018b) and Swift JO243.6+6124 (van den Eijnden 
fet al.2019). Intriguingly this latter source simultaneously launches both a wind and a jet (see Section 2.4.1). In 


addition to these studies, a larger scale search within the data of 19 ULXs, using a faster algorithm by 
(2021), indicates that rest-frame emission lines and blue-shifted absorption lines are indeed a ubiquitous feature. 

The velocity of the wind, vying, a measure of the ionisation state, € (both obtained by using consistent models for 
photoionised plasma - [Pinto et al.[2016), and an assumption for the ionising luminosity, allows the kinetic luminosity 
of the wind to be estimated from: 


Pwind = L rinna = Dabo = Vina (8) 

The filling and covering factor of the wind (Q) in the above formula remain unknown, but the high projected 
velocities (which are therefore lower limits on the true velocity) imply the kinetic luminosities for these winds are 
substantial, and may even dominate over the observed radiated power from the source (Pinto et al.|2016). It is 
worth noting that such observations are potentially in conflict with MHD models which predict lower fractions of the 
radiative power spent in launching the wind (Jiang et al./2014). Questions remain over the nature of the rest-frame 
emission lines - conceivably these may be due to collisional ionization as a result of interaction with the wind of the 
companion star (e.g. Mauerhan et al./2010) or may be associated with lower velocity winds driven 
from the outer, irradiated disc (Middleton et al.|202 1a). 

Studies have recently shifted to determining how the wind appears to change as a function of spectral hardness. 
As first reported in [Middleton et al.) (2015c), the strength of the CCD residuals in NGC 1313 X-1 appears to anti- 
correlate with increasing spectral hardness (i.e. the features are weaker when the source appears spectrally harder). 
One interpretation has been that the wind or system is precessing, or the mass accretion rate changing (the optical 
depth along our line-of-sight through the wind changing with either/both: [Poutanen et al./2007). A more recent in- 
depth study of NGC 1313 X-1 based on a series of long XMM-Newton observations, has resolved the atomic features 
in emission and absorption (using the RGS) and made a direct comparison to changes in the X-ray spectrum (Pinto| 
fet al.[2020b). These data have allowed for important progress to be made, including the discovery of a slower, more 
neutral component of the wind (moving with a projected velocity of ~ 0.08c rather than the fast wind of ~ 0.2c), 
pointing towards a complex outflow (as one might expect from simulations - e.g. ‘Takeuchi et al.]2013). These high 
energy-resolution results broadly confirm the picture from the CCD-quality residuals (as does the larger sample result 
of|Kosec et al.[2021), with the observation of an anti-correlation between spectral hardness and line strength, possibly 
driven by changes in accretion rate and/or inclination angle via precession (Middleton et al.!2015b/d). Importantly the 
correlated changes further confirm that the winds are associated with the ULX rather than some surrounding, more 
distant material (which would not respond on such timescales). The high quality data made available via long RGS 
exposures has allowed, for the first time, photoionization modelling of a ULX wind see also|Pinto| 
fet al.[2021| for a comparison to other ULXs), indicating that material in the innermost regions is potentially thermally 
unstable. In future, the launch of XRISM and Athena (with its X-ray Integral Field Unit) will allow for far deeper 
studies of ULX winds. 


2.2.2. Spectral studies of neutron star ULXs 

The discovery of several pulsating neutron star ULXs (Bachetti et al.|2014}|Israel et al./2017a\b} 
DARSE R SOTA es Table ba 
led to scrutiny of the X-ray spectrum in light of the expected differences from the standard super-critical model, as a 
consequence of the presence of the magnetosphere — accretion curtain and column — and surface. The truncation of the 
disc at the magnetospheric radius and free-fall of material onto the magnetic poles, leads to two likely configurations 
depending on the accretion rate. At very high accretion rates or low dipole field strengths, the spherization radius is 
reached at radii greater than the magnetospheric radius — if the wind is optically thick (see Vasilopoulos et al.[2019} 
Middleton et al./2019b) it will act to collimate some of the radiation within. Regardless, we expect thermal emission 
from Rspn down to Rmag, which tends increasingly towards a blackbody as Rspn tends towards Rag. At radii larger 
than Rspn, the emission is from an unmodified, approximately thin disc (thin below ~ 30% Eddington, e.g. 
fet al.[2006), although this emission is probably modified and scattered by the wind covering those regions between 
Rgpn and Rph. For strong magnetic fields or relatively low accretion rates, Rsph < Rmag and the disc emission beyond 
Rmag Originates from an Eddington—limited thin disc. In either accretion-rate regime, within Rmag, the material falls 
freely along magnetic field lines until meeting the standing shock front in the accretion column. The key characteristic 
of the magnetospheric curtain is the optical depth t which depends on the accretion rate (and therefore dipole field 


strength) and angle from the accretion disc plane. According to (Mushtukov et al.!2017) the optical depth is given by: 
TOL 39 By’ pa ) 
Ty ——— 
BO) 


where 6) ~ m — (R/Rm)!/” is the angle to the accretion column base and (0) = v(@)/c is the local dimensionless 
velocity along the magnetic field lines. This optical depth determines the shape of the thermal spectrum from the 
curtain (becoming optically thick for L39 > ay) with a range of emergent photon temperatures due to inclination 


dependence, as well as potentially diluting variability/Mushtukov et al.[2019) and the presence of cyclotron resonance 
scattering features (see the following section). A rough estimate for the temperature of the curtain is provided by 
as Tou ~ 0.5L 49! 8B,5!7m-"/4R-°!7 keV (where Re = R/10° cm and m = M/Mg) which can 
be = | keV for typical ULX luminosities, and could therefore explain one of the hard components in the spectra of 
ULXs. 
The final set of likely spectral com- 

ponents is then: at high accretion 
rates/low dipole field strength — thin disc, 
thick disc and wind, curtain and column, 
and for low accretion rates/high dipole 
field strength — thin disc, curtain and col- 
umn. Using the latter description of the A102 
inflow, |Koliopanos et al.|(2017) explored a 
the possible range of dipole magnetic g 
field strengths for 18 ULXs (containing > 
both pulsing ULXs and those without) Š% 

by 

m 


(9) 


cos 0 


under the assumption that Rsph < Rmag, 
finding values ranging from ~ 10!! — 
1013 G (although there is a tension 
between the assumption of a thin disc 10° 
and observed super-Eddington luminos- 
ity from the soft component). 

The spectrum of the pulsed emission 
can also be extracted via selecting ‘off’ Energy (keV) 
versus ‘on’ pulse times and obtaining the 
difference spectrum. The latter has been 


Figure 5 Pulse-on versus pulse-off spectra for NGC 7793 P13 from 


performed by |Brightman et al.|(2016) in , 
the case of MSD X-2, and (2018b), using data from NuSTAR. 
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in the case of NGC 7793 P13 

(Figure 5), who found the accretion column to be well-modelled by a power-law with an exponential cutoff (see 
also for a similar finding for SMC X-3). The unambiguously spectrally hard nature of the 
pulsed emission has naturally led to suggestions that other, similarly hard ULXs (e.g. Ho IX X-1: 
and IC 342 X-1{Middleton et al.[2015a) may also harbour neutron stars 
[Walton et al./2018c}|Gurpide et al.|2021). Showing whether the ‘hard excess’ is indeed due to the presence of an 


accretion column will be valuable for identifying neutron star ULXs without the need for deep pulsation searches; this 
in turn will help reveal ULX demographics and test our grasp of binary evolution. 


2.2.3. Cyclotron resonance scattering features 

In the presence of strong magnetic fields, electron and proton orbits are quantised into Landau levels; such particles 
can resonantly scatter incident photons of sufficient energy depending on the quantum mechanical cross section (see 
Schwarm et al.|2017) and leave cyclotron resonance scattering features (CRSFs) imprinted on the spectrum (although 
any subsequent scattering, for instance in an accretion curtain, can reduce their detectability - see 
(2017). The transition energy of such lines is related to the magnetic field strength and particle type according to: 


11.6 B 
tees) (ara) a a 

6.3 B 
Ee aed) free, oe (1i) 

where z is the gravitational redshift, given by: 
-1/2 
gee i (12) 
FoycC? 


where reyc is the distance from the surface of the neutron star to where the line is formed. Studies of CRSFs in Galactic 


systems (e.g. |Tsygankov et al./2006}|Jaisawal and Naik|2016) have been useful for studying the regions in which these 
lines are formed (see the review of 2019). 


reported the first CRSF in a ULX, at an energy of ~ 4.5 keV in a Chandra observation 
of M51 ULX-8 (see Figure 6); as this energy does not correspond to known instrumental features and cannot be 
readily explained as a highly blueshifted absorption line, a CRSF is a possible explanation. Given the line energy, the 
implications are either electrons orbiting in a field with a strength of ~ 5 x 10!! G (for z = 0.25) or protons orbiting 
in a field around 9x10!4 G. Given the narrow width of the line, it was concluded by (2018) that 
the feature is a proton CRSF. Combined with constraints from spectral fitting, it was concluded by 
that the dipole field is likely to be weak with the feature the result of a stronger, higher order multipolar 


(e.g. quadrupole) component closer to the neutron star (e.g. {Israel et al.[2017alb} [Brice et al.]2021), a situation now 
confirmed in the case of Swift J0O243.6+6124 

At the time of writing, there is only a single other claim of a CRSF in a ULX, in the transient system, NGC 300 
ULX-1, which was identified as a neutron star from its 31.6 second period (Carpano et al./2018), with the presence 
of a CRSF inferred from spectral modelling of the pulsed component EET the line energy 
(= 13 keV) appears consistent with the estimate of the dipole field strength from spin-up, doubts have been raised 


by |Koliopanos et al.|(2019) as to the presence of the line, as applying a spectral model based on emission from the 


accretion curtain and additional hard power-law tail, do not argue for its inclusion. 


2.2.4. Spectral evolution 

ULXs are well known to change in both brightness and in spectral shape over timescales of months or less. Most 
ULXs tend to follow a roughly predictable track with a given source becoming brighter as the hard and soft spectral 
components both increase in brightness, with the harder component tending to increase more than the soft (see e.g. 
although this depends on how the spectrum is modelled (as this delineates the hard and soft 
components). As a result, for some of the brightest ULXs, there tends to be a clear correlation between the source 
becoming spectrally harder (with more flux emerging above ~ 1 keV compared to that below) and brighter 


2018) although this is not universal (see|Gurpide et al.|2021). 
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Figure 6 Figure adapted from|Brightman et al.|(2018) showing the location of ULX-8 in M51 (left), and the CRSF 
identified at ~4.5 keV in a Chandra observation (right). 


The coupled changes in X-ray spectrum and observed luminosity can be explained by a given ULX undergoing a 
change in accretion rate (at large radii), inclination angle (due to precession) or a combination of the two 
fet al.[2015a). Under the condition that the neutron star dipole field is weak, so that Rp, < Rm, or the compact object 
is a black hole, increasing the accretion rate increases the spherization radius and the optical depth of the wind (as it 
is an integrated quantity — see{Poutanen et al.[2007) and decreases the temperature at both the spherization radius and 
outer photospheric radius (see equations 5 & 6). Accordingly, we should expect to see an anticorrelation between the 
soft component’s luminosity (which may be driven by radii interior to Rspn) and the temperature at the spherization 
radius. Identifying and interpreting a correlation (or lack of one) is not necessarily straightforward, as increasing 
Rsph, OF Moving to a more face-on orientation, leads to steadily increased beaming of the radiation interior to Rsph 
which can make unambiguous spectral modelling challenging. Conversely, the hard X-ray 
emission is not expected to be particularly sensitive to accretion rate (see Eq. |5) but is the region deepest within the 
wind cone and therefore experiences the largest amount of collimation and beaming (this is an important point when 
considering the effects on the local environment - see Section 2.6). An increase in accretion rate is only likely to affect 


the luminosity of the hard component if the wind opening angle is sensitive to accretion rate (c.f. 2014 
2019a). For a fixed inclination where sight-lines graze or intercept the wind, this picture is further complicated by 


obscuration and reprocessing (Middleton et al./2015a} although the effect of absolute collimation at any angle can be 
inferred numerically — 2017). 


Given the observed precession of S8433 (known to be highly super-critical and established to be a ULX seen edge- 
on - see Section 2.8.3), it is reasonable to expect precession to occur in other ULXs as well (for a variety of reasons, 


and with some support from observation e.g. |Pasham and Strohmayer|2013}|Middleton et al.!2015c}|/Weng and Feng 


2018). It is therefore important to consider the effect such a changing effective inclination of the super-critical flow 


might have (see |Middleton et al.|2015a/for details). Also, whilst changes in accretion rate and/or inclination lead to 


relatively predictable coupled evolution in spectrum and brightness, we also expect an impact on the timing signatures 
(see[Middleton et al.[2015a), and detected strength and speed of winds (see [Pinto et al.[2020al |2021). Such additional 
considerations are vital as they are independent of the spectral model used to describe the emission which can often 
be degenerate (see Section 2.3.4). 

The majority of the spectral evolution observed in ULXs can be fitted into the above picture. Notably the model 
can broadly explain observations of super soft ULXs, which have extremely low temperature soft components (kTpp ~ 


50 — 150eV) and very little emission above 1 keV (e.g. |Kong and Di Stefano/2003}|Urquhart and Soria|2016a 
2016) and which can be explained by an edge-on orientation and/or very high accretion rates (where the 
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wind is highly obscuring). Conversely, the brightest and spectrally hardest ULXs, tend to become slightly softer 
as they brighten, consistent with increased beaming of emission within Rspn and where 
the ULX is viewed at lower inclinations. In spite of the success of this simple model, a number of important open 
questions remain. 

One of the most obvious questions facing the unified inclination/accretion rate model of ULXs is how the presence 
of a strong dipole field affects the flow. For the most part, and as long as Ry < Rgpn, the expectation is that the 
spectrum will evolve along similar paths, as the position of the magnetospheric radius is expected to be independent 
of the mass-transfer rate for classical supercritical flows (see(Chashkina et al.[2019]for a discussion of how the location 
of Ry can be affected by advection). Rather more drastic deviations from the model above occur for those NS ULXs 
(note not necessarily ULX pulsars, as there are many reasons why a pulsation may be absent — e.g. 
which are close to spin equilibrium, as the onset of a propeller state can switch off the accretion curtain and 
column, leading to a quenching of the hardest spectral component (while the other components remain ‘on’). There 
are already hints of propeller states occurring in one ULX pulsar (Tsygankov et al.|/2016b) and such states might 
account for some of the more highly variable ULXs (e.g. although a 
propeller state is ruled out for recent, large-scale changes in X-ray flux observed in NGC 300 ULX: 
and NGC 7793 P13: |Fuerst et al.}2021). 

One of the most intriguing observations is that the shape and normalisation of the emission above 10 keV tends 


to stay remarkably constant in a number of bright ULXs as their spectra evolve (see 2020). This can be 
explained so long as the innermost regions are maximally beamed and any changes in accretion rate do not change the 


opening angle of the wind cone appreciably (although see Jiang et al.[2019alfor indications that the cone’s opening 
angle may well respond to larger changes in accretion rate). An alternative possibility is that the emission at these 
energies originates from a stable accretion column, and applying a phenomenological model for the column can indeed 
provide a good description of the hard emission [Walton et al.[2018c). 

There are important questions remaining over how the individual spectral components fit into the above picture, 
typically explored through use of the temperature-luminosity plane. Such analyses have a long history in the study 
of ULXs, with indications of a strong negative L-T correlation (L œ T~>°) reported across a sample of ULXs (e.g. 
[Kajava and Poutanen]2008] |2009). Given the expectation that the ULX population is heterogeneous, it is important to 
also study sources individually; this is becoming increasingly feasible with intense observing campaigns. 

examined nine observations of NGC 1313 X-1 (combining XMM-Newton with NuSTAR data) 
and fitted the spectra with a two-component thermal model of a disc black-body at low energies and an advection 
dominated disc at high energies (and an additional component to account for the hard excess). From this they found a 
positive (and relatively steep) correlation between the temperature and luminosity of the hotter, advection dominated 
disc (the soft component has a flatter, possibly negative correlation). This correlation is further subdivided into high 
and low luminosity branches and may provide an intriguing test of what processes drive the changing appearance of 
ULXs. A similar analysis has been performed by|Robba et al.|(2021) for the ULX pulsar NGC 1313 X-2, although a 
second luminosity-temperature branch is not as evident as for NGC 1313 X-1. As shown in several works 
[Gúrpide et al./2021), the presence and nature of any correlation (at low or high temperatures) can be 
affected by the choice of spectral model, reinforcing the important role other methods can play in understanding the 
nature of the flow and the components in the X-ray spectrum (see Section 2.3.4). 


2.2.5. Open questions and future directions 
Important questions remain regarding the modelling of ULX spectra, including: 


e the apparent stability and similarity at energies above the spectral break (Walton et al./2020) 


e the presence of fainter spectral components which may correspond to bremsstrahlung emission, as seen from 


SS433’s jets (e.g. [Walton et al./2015) 


e the trends of luminosity versus temperature for the individually modelled components (Kajava and Poutanen 
2009} [Walton et al.|2020). 


e the impact of spin and advection on the hard spectral emission 


e the prevalence (or lack-thereof) of CRSFs 


Many of these issues may be resolved as self-consistent spectral models derived directly from 3D simulations of 


black hole and neutron star ULXs become available, (e.g. |Ohsuga and Mineshige}201 1} 2014} [Sadowski 
et al.|2015}|Sadowski and Narayan|2016}|Narayan et al.|2017||Takahashi et al.|2018}|Kawashima and Ohsuga|2020). 


2.3. Timing analysis 


X-ray spectroscopy has allowed for a broad categorisation of ULXs, but as with other accreting systems, the timing 
properties (over short and long timescales) are just as important for testing models and for revealing new behaviours. 


2.3.1. Short timescales 

Because of its large effective area, XMM-Newton has been the favoured instrument for studying the short (intra— 
observation) timescale variability of ULXs. The simplest analysis uses the Fast Fourier Transform (FFT) of the 
lightcurve to obtain the periodogram (or, if averaged, the power spectral density — PSD). As this is typically normalised 
to be in (rms/mean)?/Hz units, the integral of the PSD gives the fractional root mean square (rms), which, when 


Poisson noise corrected and plotted against energy, yields the rms spectrum (see|Uttley et al./2014|for a comprehensive 
review). 


The PSDs extracted from single observations of ULXs are generally described by either a red—noise power law 
(power scaling as v-! — v~?) or a broken power-law (e.g. 
Atapin et al.[2019), sometimes designated as the broadband noise. Binning 
the rms reveals a tentative linear relationship with X-ray flux 
Garcia et al.|2015) as expected from multiplicative propagation of fluctuations in accretion rate (Lyubarskii 

[Arévalo and Uttley|/2006). observed that the variability of some ULXs 
appeared very low, while conversely others appeared extremely high (up to 30% fractional rms - [Middleton et al.] 
[2015a). This can be explained via either the dilution of intrinsic variability from a neutron-star accretion column as 
photons scatter inside the accretion curtain before eventually escaping (Mushtukov et al.2019) or by mass loss in 
the wind, which may remove variability from the inflow (Middleton et al./2015a). The substantial amounts of short 
timescale (< 10 ks) variability in those spectrally softer sources — in which there is strong evidence for mass-loaded 
winds (e.g. NGC 5408 X-1 and NGC 6946 X-1) — could arise from stochastic obscuration by clumps in the wind 
Sutton et al.|2013), produced via radiative hydrodynamic (e.g. or Kelvin 
Helmholtz instabilities. 

In addition to the broadband noise, there have been several claims of quasi-periodic oscillations (QPOs — with 
power concentrated over a small frequency range) at 10s of mHz (NGC 5408 X-1: 
NGC 6946 X1 NGC 1313 XE Pa s ATOS MED XT 
IC 342 X-1: as well as at higher (~Hz) frequencies in 
M82 X-1 (Pasham et al.2014). There have been various attempts to use QPOs to obtain the accretor mass, based 


upon a simple frequency scaling (e.g. {Strohmayer and Mushotzky||2009 2015), which tend to give 


masses typical of IMBHs. However, while QPOs are detected from super-critically accreting systems such as TDEs 
(e.g. and in GRS 1915+105 (the well known 67-68 Hz QPO: [Belloni and Altamirano[2013), it is 
not clear that there are well-studied objects with well-understood signals available for direct comparison. At present, 
the origin of ULX QPOs remains unclear, although models invoking modulations of accretion rate 
or precession of the radiation-pressure supported disc — which can also explain the lag seen around the QPO 
frequency (De Marco et al.[2013]— see Section 2.3.4) have both been put forward (Middleton et al.[2019b). With the 
high throughput of forthcoming instruments (namely ESA’s Athena), studies of ULX QPOs will become considerably 
more revealing. 


2.3.2. Pulsations 
Following the discovery of an X-ray pulsation from M82 X-2 (Bachetti et al./2014), one of the most important 


tools now regularly applied to ULX data is the accelerated pulsation search (see e.g. 2001 
2002} [Andersen and Ransom|2018). This can account for the orbit of the neutron star about the system barycenter, 


which acts to smear the power in the pulsations over neighbouring Fourier frequency bins. Accelerated searches 
within both new and archival data have led to the discovery of several more pulsating ULXs (NGC 5907 ULX-1: 


Israel et al.|2017a; NGC 7793 P13: [Fürst et al.]2016} [Israel et al./2017b| SMC X-3:|Tsygankov et al.|2017} NGC300 
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Figure 7 From (2020) showing the positive glitch (with a change in frequency Av ~ 8x107!!Hz/s) in 


M82 X-2, revealed by extensive monitoring with NuSTAR. 


ULX-1: |Carpano et al./2018)— first identified as a supernova by NGC 1313 X-2: 
M51 ULX-7: RX J0209.6-7427: 
2020) and Swift JO243.6+6124: |Kennea et al.|2017). 

The observed pulsing ULXs typically have periods ranging from ~ 1 - 30s (see Table B). As is apparent from 
Tables |2| the secular P in these systems (i.e. AP/AT where AT is typically months or longer — except in the case 
of NGC 300 ULX-1) reveals a collection of objects which tend to spin up at extremely high rates (notably NGC 


300 ULX: 2018) when compared to Galactic X-ray pulsars (although M82 X-2 is observed to spin 


down between 2014 and 2016: 2020). The secular spin up/down rates may differ greatly from the 
‘instantaneous’ P obtained over the course of an observation due to the influence of orbital dynamics and changes 


in accretion rate between observations (Israel et al./2017a 2020). In addition to these long timescale 
trends, glitches (perhaps associated with the coupling of differentially rotating inner crust — and possibly core — to 


outer crust: [Anderson and Itoh]1975) have now been observed in M82 X-2 (Bachetti et al.|2020|— Figure 7). These 


imply a change in spin frequency v/v ~ 1075, similar in magnitude to glitches seen in magnetars and a small number 
of acereting neutron stars 
et al.[2017). Interestingly, there are two (or possibly three) anti-glitches in NGC 300 ULX-1 (Figure 8, 
2019), where the neutron star undergoes a dramatic spin-down dv/v > —10-+. This may be the result of the enormous 
spin-up rate creating a sizeable lag between the rotation of the crust and superfluid interior (Ray et al./2019) and is 
supported by the lack of corresponding radiative changes during the glitch. Only spin-up glitches have been observed 
in non-accreting pulsars, but accreting pulsars and magnetars show both spin-up and spin-down 


glitches (e.g. Dib and Kaspi/2014). The anti-glitches seen in NGC 300 ULX-1 are the largest yet-observed in any 
kind of pulsar system. 


The known ULX pulsars have pulse fractions (defined as [Pulse max - Pulse min]/[Pulse max + Pulse min]) which 
increase with photon energy (as seen in Galactic X-ray pulsars: Ea and take values 
from <5-100 % depending on the observation, source and energy band (see Table|2|for more details). Notably, Figure 
9 shows that several of the ULX pulsars have an a highly sinusoidal pulse profile (although both NGC 300 ULX- 
1 and SMC X-3 are notable exceptions), unlike that expected from a pencil beam geometry. This suggests a large 


emitting region, either a fan-beam geometry (Basko and Sunyaev| 1976) or accretion curtain, the latter reprocessing 
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Figure 8 Data from (2019), highlighting the glitches in the ULX pulsar NGC 300 ULX-1. The upper 
(orange) points show offset residuals from the best-fitting model of the frequency evolution without glitches. These 


anti-glitches (with ôv/v > —1074) are the largest yet observed in any pulsar system, and may be a result of the lag 
between the rotation of the crust and superfluid interior, driven by the enormous spin-up rate (Table[2). 


(and therefore potentially diluting) the pulsed emission (Mushtukov et al./2017||2019). 


We can obtain a rough indication of the statistical significance of a given pulsed signal (with a width Anu, fractional 
rms, r — related to the pulse fraction — in an observation of length T, from a source of count rate 7, with background 
count rate B) from no = 1/2 UP /(I + B)]r? VT/Av noting that a true measure of the significance 
relies on considering both the noise and the number of independent free trials). Of the ~ 1800 known ULXs (e.g. 


2022), the vast majority do not have have sufficient data quality for deep pulse searches (down to small 
r). In those with sufficient data quality and not already identified as PULXs, pulsations appear to be absent down to 


the detection threshold (e.g. |Doroshenko et al./2015 2017b). 
There are several reasons why pulsations may be absent in a given ULX or in a given observation. Pulsations 
are observed to be transient in the known ULX pulsars (and indeed in Galactic accreting X-ray pulsars as well, e.g. 


2000) for reasons not yet well understood (see|Bachetti et al.|2014 2017a\|Sathyaprakash 
et al.|2019b}|/Rodriguez Castillo et al.|2020 2020), and conditions during a given observation may not be 


favourable for detecting the pulsed emission from the column (and/or curtain). As an example, we may be prevented 
from observing the pulsations by dilution by a stronger stable component in the spectrum (Walton et al./2018b} e.g. 


if the disc/wind has precessed or accretion rate has changed|Pasham and Strohmayer|2013}|Middleton et al.|2015a 


2018|/2019b). This can also happen if the neutron star has free-body precessed (Fürst et al.|2017) so that the pulses are 
considerably weaker or beamed away from us. It may also be that we will never observe pulsations in some neutron 


star ULXs as we may never observe at suitable angles to view the accretion column. Alternatively the field may have 
been suppressed by accretion (e.g. [Urpin and Geppert|1995), the neutron star and inner disc may be aligned (King and] 
Lasota|2020), scattering in the wind cone may have suppressed pulsations on timescales shorter than the light-crossing 


time (Mushtukov et al./2020), or some ULXs may simply harbour black holes (see|Middleton and King}2017| Binary 
evolution would seem to make the latter inevitable to some degree — (Wiktorowicz et al.|/2019). 
2.3.3. Long timescales 


It is well-accepted that LMXBs undergoing recurrent outbursts (see|Hameury and Lasotal2020) may appear as 
ULXs for relatively short periods of time (e.g. |Middleton et al./2013 2018} |van den Eijnden et al. 
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Figure 9 Pulse profiles for the PULXs shown in Table |2| extracted from data taken by XMM-Newton (NGC 5907 
ULX-1, M51 ULX-7, NGC 7793 P13, NGC 300 ULX-1, NGC 1313 X-1) and NuSTAR (M82 X-2 and SMC X-3). 
Values for the count rates are not provided due to bandpass and detector differences, however, errors are provided for 
the purposes of comparing the level of constraints within each phase bin. 


2018). Many bright ULXs are persistent (although we note the discovery of a new transient ULX in M51 reaching 


10*° erg/s: [Brightman et al.[2020a), but their luminosities can fluctuate by orders of magnitude, periodically/quasi- 
periodically (on timescales of 10s of days: or seemingly 
at random. There are several possible explanations for such variations: changes in accretion rate, driven by effects in 
the outer disc (see Fig. [T9] Sec. and[Hameury and Lasota|2020} [Middleton et al.[2021a), an accreting neutron 
star entering a propeller phase, orbital variability as a neutron star passes through the decretion disc of a Be star, and 
precession of the disc/wind under external torques. In a number of the pulsating ULXs, the orbital period can be 
deduced from the Doppler modulation of the pulsations; this gives periods of a few days (see Table B. with NGC 
7793 P13 having a proposed period of 64 days (Motch et al./2014]although see [Fuerst et al-[202i]for details regarding 
the periodicities in this ULX). Other than NGC 7793 P13, the ~ day orbital periods are considerably shorter than the 
“super-orbital’ periods typically revealed by Swift monitoring but in some cases observed in the optical/UV (e.g. NGC 
7793 P13: Fuerst et al.|2021). Both sets of periods are reported in Table[2| 

Precession has been proposed as one means by which to explain the ~ 10s of days super-orbital variability of 
ULXs by a number of authors. 2013) suggested that the changes in the spectrum of M82 
X-1 associated with a 62 day period (Kaaret et al.|/2006} |Kaaret and Feng|2007) could result from changes in the 
projected area of the inner disc (perhaps driven by a radiative warp: ; (2015c) suggested 
that changes in the strength of atomic features, now known to be associated with the wind (Pinto et al./2016 [2020b), 
could result from precession of the disc and wind in NGC 1313 X-1, and{Luangtip et al.|(2016) appeal to a combination 
of precession and changes in accretion rate in order to explain the changing appearance of Ho IX X-1. 

There are a number of possible mechanisms for driving a changing view of the inner regions of ULXs, including 
precession of a neutron star’s magnetic dipole field (Lipunov and Shakura] 1980): 
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where u30 = BnR, I45 is the neutron star moment-of-inertia, ¢ is the angle between the normal to the disc plane 
and the spin axis, and ¢ is the angle between the magnetic axis and the spin axis. To match observations of ULX 


superorbital periods of months or less, field strengths > 10'4G are required 
etal. 2020a). 

Should the accretion flow be misaligned with the neutron star or black hole spin axis (which seems quite plausible: 
then frame-dragging of the fluid induces a torque. Similar to the hot inner flow of X-ray 
binaries, a thick super-critical disc/wind is expected to have an effective viscosity parameter a < H/R. The frame- 
dragging induced torque is then expected to be communicated as a bending wave from the inner edge of the flow out 


to the photospheric radius (Middleton et al./2018}/2019b) which may then induce solid body precession (Fragile et al. 
2007) on a estimated timescale of: 
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for the wind (here r is in units of gravitational radius - GM/c? — and a, is the dimensionless spin parameter). Both 
periods are sensitive to the mass accretion rate. This can be related to properties of the X-ray spectrum and may tell us 
the nature of the accretor (neutron star or black hole, [Middleton et al./2019b) without the requirement for pulsations 


which are themselves transient or may be diluted beyond detection (Bachetti et al.|2014}|Mushtukov et al./2020). 


If the neutron star is asymmetric about its rotation axis, free-body precession may also occur on a timescale of 


(Jones and Andersson|2001): 
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where y is the angle of the misaligned distortion of the neutron star relative to its rotation axis, and € is the measure 
of distortion, related to the surface magnetic field strength through (Lander and Jones|2009): 
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which, notably, [Vasilopoulos et al.|(2020a) used to estimate the dipole field strength in M51 ULX-7. 

Besides the above torques, other effects may warp and precess the disc, e.g tidal torques from the secondary 
star (although these are unlikely to be strong in cases with known ~ day orbital periods, e.g. |Bachetti et al./2014), 
magnetic warping Pfeiffer and Laif2004), and radiative warping (Pringle]1996). In the latter case, where 
super-critical flows — specifically a large-scale height disc and optically thick wind — are present, it is not clear how 
efficiently the outer disc is irradiated (see also [Middleton et al./2021a). But if the dipole field is strong and/or the 
accretion rate low, irradiation may well be possible, and similar warps and precession to that seen in Her X-1 are then 
quite plausible (see [Petterson]1977| 2012). Regardless of the mechanism driving the precession, 
it is possible to build a simplified model of the motion of the cone and then numerically extract the X-ray lightcurve. 
(2017) developed a model which can produce a variety of lightcurve profiles (available within xsPEc: 
as ULXLC), under the assumption that the X-ray source lies at the apex of a cone of outflowing, 
totally opaque, achromatically scattering plasma. Applying this model to the lightcurve of NGC 5907 ULX-1 implies 
a small total opening angle (~ 10 — 15 deg) and substantial beaming (amplification by a factor 60-90) although it is as 
yet unclear whether this is in conflict with scattering reducing the observed pulse fraction (Mushtukov et al.[2020). 

It is interesting to consider how the different physical origins for superorbital periods can be distinguished. A 
natural method would be to consider how the timescales respond to the accretion rate and spin-up/down of the neutron 
star. Equations 14 & 15 are highly sensitive to the accretion rate (far more so than the spin) and so, where the accretion 
rate is variable (on long timescales), we would expect variability in the superorbital period (with correlated changes 
in the X-ray spectrum). Conversely, magnetic and free-body precession is more sensitive to the spin than accretion 
rate, and the superorbital period should therefore track changes in the neutron star’s spin over time. 

It is also worth noting that plotting the superorbital period versus the orbital period for Galactic systems and 
those ULXs with known periods (see Table [2| implies a potential correlation (see [Townsend and Charles|{2020), 
although care must be taken in interpreting any correlation between logarithmically scaled values which include 
systems known to be of very different types. It is plausible that the superorbital periods are in some cases responding 
to tidal interactions, while the mass supply to the compact object may correlate with the binary period, leading to a 
positive correlation through e.g. equation 15. This currently remains an open question but we can expect progress as 
the detected number of constrained superorbital periods increases. 

In addition to the smooth, periodic/quasi-periodic variations seen in the X-ray and optical bands, there are a variety 
of changes in the lightcurves of ULXs which can help reveal the nature of the system. Dips are found in NGC 55 
ULX-1 (Stobbart et al./2004); NGC 5408 X-1 (Grisé et al.|2013), NGC 247 ULX-1 (Alston et al./2021), M51 ULX-7 
(a known ULX pulsar, {Vasilopoulos et al.[2021), and J125048.6+410743 (Lin et al.J2013). Significantly, these ULXs 
tend to have relatively soft spectra (M51 ULX-7 is somewhat harder), possibly implying a correlation with sight-lines 


intercepting regions in the wind subject to instabilities and clumping (e.g. |Middleton et al.2011||2015a 
et al./2013). 


Whilst dips are seen in the lightcurves of a small number of ULXs, eclipses seem to be mostly absent from the 


currently known ULX population (Middleton and King|2016). This would make sense if edge-on ULXs emit around 
the Eddington limit (Poutanen et al./2007) implying sub-ULX luminosities in many cases (Zhou et al./2019). However, 
there have been reports of eclipses in two ULXs in M51 (Urquhart and Soria/2016b; 2018). These give 


constraints on the angle between the companion star’s orbit and the X-ray bright inner regions (noting that these do 
not have to be aligned in a single plane). If there is an estimate of the mass of the donor star this then constrains the 


compact object mass (Wang et al.|2018). 
Neutron stars probably account for the majority of observed ULXs (King and Lasota)2016} |Middleton and King 
2017 2017} 2018c 2019). Such systems are likely to approach spin 
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equilibrium between gravity and centrifugal forces during their lifetimes, where the co-rotation radius approaches the 
magnetospheric radius. Whilst the effective torque vanishes when the two are close (Dai and Li] /2006), should some 
process push the magnetospheric radius outside the corotation radius, a propeller state ensues where the emission 
from the accretion column and curtain — but importantly not from radii down to the magnetospheric radius — is 
switched off. In systems where the magnetospheric radius is large, the result would be a severe drop in flux (although 
centrifugally driven winds from around the magnetospheric radius would remain). Such propeller states, widely 
observed in Galactic systems (e.g. have been invoked to 
explain the large changes in luminosity seen in M82 X-2 (Tsygankov et al. . However, in NGC 7793 P13 
(Fuerst et al.|2021) and NGC 300 ULX-1 (Vasilopoulos et al./2019), the spin has continued increasing at the same rate 
when the source is considerably fainter in the X-rays, implying that some other mechanism — possibly obscuration — 
drives the change, rather than a propeller state. (In M51 ULX-7 it is at yet unclear whether a similar X-ray off state 


implies a propeller state or is associated with precession: | Vasilopoulos et al.|2021). The predicted large-scale changes 


in X-ray flux in the propeller state have motivated searches for pulsar candidates among known ULXs. 
searched within those ULXs in the 3XMM DR4 catalogue finding five objects which 
underwent a factor of 10 or more change in flux (with one propeller candidate: M51 ULX-4). (2020) 
compiled XMM-Newton, Swift and Chandra flux lightcurves of 296 ULXs from the (2019) catalogue, 
and found that 25 out of 296 ULXs changed flux by more than a factor 10 (based on the available sampling), 17 of 
which show a bi-modal flux distribution, potentially consistent with sources at spin equilibrium. 


2.3.4. Spectral-Timing studies 

As well as standalone studies focussing on either the spectral evolution or timing characteristics of ULXs, there 
are an increasing number which bring the two together. Such a spectral-timing analysis can be performed in the 
time or Fourier domain in a given energy band (e.g. the fractional rms or rms spectrum: 
Gilfanov et al. and the covariance spectrum: [Wilkinson and Uttley/2009), or via lags between bands (via the 
CCF or Fourier lag spectrum) as a function of energy or frequency (see for a comprehensive review). 
Spectral-timing studies require high signal-to-noise or large amounts of variability (such that the Poisson noise is not 
dominant), making such analyses challenging for ULXs. However, there have been a number of successful attempts 
which have provided otherwise inaccessible insights. 


Covariance spectra have an advantage over rms spectra; as the Poisson noise between 
bands is uncorrelated, when the variability is weighted by the linear coherence between bands (see 
[Uttley et al.[2014), the constraints on the variability in each energy bin are improved. 
(20 15a) presented the covariance spectra for a sample of variable ULXs with relatively high data quality and showed 
that the variability on relatively short (S 1 ks) timescales matched well to the shape of the hard component (seen 
also the studies of the covariance spectra of NGC 55 ULX-1: and NGC 1313 X-1: 
(2020). This can be readily explained if the variability emerges from only this component on the timescales being 
studied. However, it is less clear why spectrally-harder ULXs (where this component is dominant) appear to be less 
variable on similar timescales (e.g. |Heil et al.]2009). One possible explanation is that the variability is extrinsic and 
being driven by obscuration by optically thick clumps of wind (Middleton et al./2011), probably created due to fluid 
instabilities (Takeuchi et al.|2013). Clumps may also then be responsible for dips seen in the softer ULXs (e.g. NGC 
55 ULX-1: |Stobbart et al.|2004/ and eet pate ale which in the standard paradigm are viewed at 
higher inclinations (Poutanen et al.|2007}|/Middleton et al./2015a). 

The covariance spectrum can also be used to isolate the non-variable component in ULXs. This can be modelled 
and better compared to theories for emission at soft energies (typically either a thin disc beyond the magnetospheric 
radius or emission from close to the spherisation radius). This can be complex, depending on the role of variable 
absorption, but it is at least possible to compare luminosities, and rule out the presence of a thin disc. This was 
demonstrated in [Middleton et al.| (2019a) (Figure 10), and allowed for independent constraints to be placed on the 
dipole field strength in M51 ULX-7 (see also discussions in and|Christodoulou et al.]2019). 

Time (or phase) lags between bands may be extracted in the time domain via the CCF (as shown to great effect 
when studying a variety of features in XRBs, e.g. lags corresponding to propagation in jets: or via 


the Fourier cross-spectrum (see|Uttley et al./2014/|for a detailed step-by-step guide). The lags are related to the nature 
and geometry of the accretion flow via the transfer function and are imprinted by propagation, and the light-travel time 


through reverberation or transmission (the latter affected by the amount and state of intervening gas). Lags in ULXs 
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Figure 10 From|Middleton et al.|(2019a) showing (left-hand panel) the covariance spectrum of M51 ULX-8 (in green) 


with the time-averaged Chandra data (in black). The covariance picks out the linearly correlated variability (see 


Wilkinson and Uttley}2009) and allows the spectrum to be decomposed. In this case, the additional soft component 


(shown in red in the middle and right-hand panels) is above the Eddington limit for a neutron star and so inconsistent 


with a thin disc. This allowed (Middleton et al.) (2019a) to place indirect constraints on the dipole magnetic field 


strength (a higher-order multipole field is not excluded). 


were first identified in a 2006 observation of NGC 5408 X-1 (Heil and Vaughan|2010) between a soft (0.3 - 1 keV) and 
hard (1 - 10 keV) band. This ‘soft lag’ is ~ a few seconds in at frequencies of ~ 
10 mHz, is ~x 150s at 0.1 mHz in NGC 1313 X-1 (Kara et al.|2020), and is a few hundred seconds at around the same 
frequency in NGC 55 ULX-1 (Pinto et al./2017). Soft lags have been interpreted as reprocessing in distant material 
(presumably the wind) and intriguingly appear around the QPO frequencies seen in NGC 5408 X-1 
(2013). Should the lag be imprinted as a consequence of resonant line scattering or Compton down-scattering, the 
magnitude can provide an estimate of the optical depth and thickness of the wind. An alternative explanation is that 
hard photons created in the accretion column are Comptonised by a surrounding envelope (Mushtukov et al./2017) 
and a lag is imprinted due to the light travel time across the curtain and the time taken to CS e 
let al.[2019). This picture would imply very large or dense curtains in the case of ULXs such as NGC 55 ULX-1. 


2.4. Multi-wavelength studies 


ULXs are defined (albeit empirically) by their X-ray luminosity, although much of the gravitational energy dis- 
sipated in the accretion flow may emerge in other wavebands through reprocessing or as kinetic luminosity in the 
wind and/or jet. It is therefore crucial to obtain as broad an SED as possible, and below we discuss insights obtained 
in several key bands (note that we discuss observations of bubble nebulae surrounding ULXs separately in Section 
2.6.1). 
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2.4.1. Radio 

Radio monitoring campaigns and cross-matching of X-ray catalogues with radio surveys 
(Pérez-Ramirez et al./2011) have returned limited numbers of ULX radio counterparts, perhaps suggesting that any 
jets are inherently weak or transient in nature. Targeted observations have tended to focus on the brightest ULXs 
(HLXs at ~ 10+! erg/s) as these remain the most promising candidates for hosting IMBHs other than AGN in dwarf 
galaxies (Farrell et al.[2009). Their suggested sub-Eddington accretion rates would naturally lead to the expectation 
of flat spectrum radio emission associated with a compact jet (if in an analogous hard X-ray state as seen in X-ray 
binaries and AGN). Placing such HLXs on the radio—X-ray fundamental plane 
opens up the possibility of constraining the compact object mass, and this approach has 
been adopted for a number of ULXs (NGC 4088-X1: NGC2276-3c: 
N5457-X9: IC 342 X-1: 


Confirming the presence of a jet in ULXs can be more difficult given the presence of a radio bright nebula in 
some cases. These can be up to 100s of pc across: 2013b), potentially analogous to the W50 nebula 
surrounding SS433 (e.g. |Pakull et al.|2006}|Pakull and Grisé|2008} [Russell et al.|2011)|Cseh et al.|2012), and energised 


by termination shocks between the surrounding gas and outflows. In only a single bright ULX — Ho II X-1 — has the 
nebula emission been resolved by using VLBI into discrete core and lobe emission although 
there are also twin lobes either side of a bright radio core in the IMBH candidate NGC 2276-3c, Ser 
(2015). The fading of Ho II X-1’s lobes over a period of around 2 years indicates that the jets are sporadic and may have 
more in common with the discrete ejections in near-Eddington LMXBs, placing a limit on the mass of the compact 
object in Ho II X-1 of > 25 Mo {Cseh et al.J2015). 

Jet-inflated nebulae may provide one of the few ways to locate edge-on ULXs. A key example is the radio-optical 
nebula in M83 (30 times more luminous than W50 at 5 GHz) and collimated jets in NGC 7793 (S26: 
Soria et al.|2010). In both cases, the estimated jet power lies above 10*° erg/s, but the observed X-ray luminosity 
to-date is 3-4 orders of magnitude lower, and the sources would therefore not qualify as ULXs. These systems are in 
many ways similar to SS433 (albeit with even higher jet kinetic luminosities) and provide ample evidence for similar 
edge-on ULXs being located in nearby galaxies. 

As well as jets associated with highly super-Eddington rates of mass transfer, transient jet ejections have been 
observed in some Galactic XRBs moving to less extreme, near-Eddington accretion rates (Fender et al.|2009). These 
have been successfully detected in nearby galaxies and can be associated with the appearance of a low luminosity 
ULX (e.g. M31: {Middleton et al.[2013). It is plausible that the transient jets detected from HLX-1 (Webb et al.[2012) 
suggest the source was moving to near-Eddington rates, which, at a luminosity of > 10’? erg/s would imply a BH mass 
of ~10,000 Mo, consistent with the mass estimate from modelling the X-ray spectrum with relativistic disc models 
(Davis et al.|2011). In future, with the advent of sensitive all-sky radio monitoring (SKA), the detection of such radio 
events should become commonplace. 

Finally, it is intriguing to note that a jet event has also been detected from the transient X-ray pulsar, Swift 
J0243.6+6124 — which contains a neutron star with B > 10! G (Doroshenko et al.|2018) — when accreting above its 
Eddington limit (van den Eijnden et al./2018). At peak radio brightness, the flux density from the jet was < 100uJy 


at 6 GHz (at an assumed distance ~ 7 kpc —|van den Eijnden et al./2018). Similar jets are therefore unlikely to be 


discovered in nearby galaxies (although there is the theoretical potential to drive more powerful jets through tapping 


the neutron star’s angular momentum: |Parfrey et al.|2016). 


2.5. Infra-red to ultraviolet 


Generally speaking, observations of accreting systems at energies intermediate to those of X-rays and radio allow 
us to study the impact of reprocessing, the highest energy populations of electrons in jets and the nature of the 
companion star. In ULXs, such studies are of vital importance for constraining the true energy budget and age of the 
system, for understanding the impact on the surroundings and how such high rates of mass transfer actually occur. 


2.5.1. Infrared 
Campaigns studying ULXs in the IR have focused mostly on the nature of the donor star, using a combination of 


imaging and spectroscopy from ground-based instruments. Initial imaging studies by (2014) identified 
a number of possible red supergiant (RSG) companion stars (11/62 ULXs in 37 galaxies within 10 Mpc) from a 
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combination of data from WHT/LIRIS, MMT/SWIRC and VLT/ISAAC. later performed NIR 
spectroscopy using Keck/MOSFIRE (Multi-Object Spectrometer for Infra-Red Exploration), identifying candidate 
counterparts for two ULXs (one in NGC 925 and one in NGC 4136) with RSGs. Building on the work of 
(2016), a large study by [López et al.[2017||[2020|made a systematic search for NIR counterparts to 113 ULXs, 
through use of a combination of instruments (imaging via LIRIS/WHT, SOFI/NTT and WIRC on the Palomar Hale 
5-m telescope, and spectroscopy with Keck/MOSFIRE). This study yielded candidate counterparts for 38 ULXs, 5 
of which are identified as RSGs. The numbers obtained from ground-based studies match well to mid-IR studies 
performed using Spitzer, where|Lau et al.|(2019) located 12 candidate RSGs in a sample of 96 ULXs within 10 Mpc. 
These results point towards a highly probable association of RSGs with some ULXs, and an additional — somewhat 
less considered — form of mass transfer, i.e. that of wind-fed Roche-lobe overflow (e.g. 
Mellah et al.]2019) which has now been incorporated into models for binary population synthesis (Wiktorowicz et al.| 
2021). 

There have been a number of focused campaigns to explore the brightest ULXs, and specifically the PULXs (de- 
termining the companion star in the latter being of obvious value when considering the evolution of these systems). 
Unfortunately for most PULXs the IR fields are extremely crowded, precluding the identification of a unique coun- 


terpart (Heida et al.|2019a). However the companion of the SN imposter, NGC 300 ULX-1 has been identified as an 
RSG (Heida et al.|2019b). 


As well as emission from the companion star, it is quite possible that synchrotron emission from jets and from 
circumbinary dust could both contribute to observed flux in the IR band. (2019) identify a mid-IR excess and 
a probable circumbinary disc in several ULXs (including Ho II X-1), and associate variable mid-IR emission from the 


X-ray spectrally hard ULX, Ho IX X-1, with a jet (as previously suggested by |Dudik et al.!2016)and|Sathyaprakash 


et al.|2022)in the case of the pulsating ULX, NGC 1313 X-2). 
Indirect constraints on the SED which photoionizes gas surrounding the ULX (see Section 2.6) can be extracted 


through modelling IR lines. (2010a) identified Ho II X-1’s nebula (thought to be inflated via jets) 
in Spitzer observations and, in (2010b) later performed photoionization modelling to suggest that the 
lines are consistent with irradiation by the observed soft X-ray and UV emission, requiring little to no beaming (of 
these components). This is completely consistent with the picture where the soft emission is less beamed than the 
hard X-rays. The MF16 nebula surrounding NGC 6946 X-1 has also been detected by Spitzer, with photo-ionisation 
modelling by|Berghea and Dudik|(2012) indicating the possible presence of an O-type supergiant donor star. 

Roche lobe, this mass transfer mechanism known as ’wind Roche lobe overflow” can remain stable even for large 
donor-star-to-accretor mass ratios 


2.5.2. Optical 

Identifying a unique optical counterpart to a ULX potentially gives a powerful probe of the system age via binary 
evolution arguments. If atomic lines can be isolated and associated with the companion, then subsequent studies 
may allow the mass of the compact object itself to be constrained by standard radial velocity techniques. There are 
a number of issues which make such investigations challenging. Given the distance to many ULXs, there may be 
several optical counterparts associated with the X-ray astrometric position (e.g. (Gladstone et al./2013). The distance 
can make lower mass companion stars undetectable, perhaps misleadingly favouring an association with any high 
mass star within the neighbourhood of the ULX. Further, the photometric signal of the star may be contaminated by 


direct emission from the wind photosphere depending on the accretion rate (Poutanen et al.|2007), and emission from 
the uncovered thin disc which may be increased by irradiation (should the geometry, advection and wind optical depth 


permit this —[Sutton et al.[2014} [Yao and Feng]2019). Finally the spectroscopic signatures (e.g. He II emission lines) 
may also be contaminated or dominated by the wind (Fabrika et al.|2015). 

Despite these challenges, there have been a number of attempts to find and study the optical counterparts of ULXs 
(besides the results discussed in the IR band above). Many studies (too many to recount in this review) over the last 20 
years have targeted individual ULXs (see e.g. [Roberts et al.[2001), whilst others have aimed at describing the larger 
population. As important examples of the latter, the sensitivity and angular resolution of HST allowed|Roberts et al.| 


(2008) to locate counterparts to four out of six ULXs, whilst {Tao et al.](2011) studied the optical counterparts of 13 


ULXs, finding the optical flux and colour to vary on the timescale of days to years, ruling out changes in extinction as 


the origin. (Gladstone et al.|(2013) utilised HST photometry to study 33 ULXs with Chandra positions, locating 13 + 
5 counterparts and finding O—star companions to be rare. 


23 


Where unique counterparts have been identified, there have been attempts to constrain the mass and species of 
compact object via dynamical arguments (most using the broadened He II 14686 emission line): 


e Using Gemini observations, (2011) attempted to constrain the mass of the compact object power- 
ing NGC 1313 X-2 (which is now known to be a neutron star: |Sathyaprakash et al./2019b) but were unable to 


do so due to the absence of sinusoidal motion in the line velocities 


° (2013) were unable to locate an orbital period, but placed a limit on the mass (< 510Mọ) of the 
compact object in NGC 5408 X-1 from a limit on the semi-amplitude in the absence of periodic modulations, 


and under the assumption that the He II line is produced in the irradiated outer accretion disc. 


e (2013) used the He II 24686 line to place a lower limit on the mass of the compact object in M101 
ULX-1 of 5Mo. 


e (2014) placed an upper limit on the mass of NGC 7793 P13 (< 15 Mo) by modelling the ~ 64 d 
period in the V and U band lightcurves as irradiation of the well-identified B9Ia companion star; this result is of 
course consistent with the discovery that a neutron star powers this ULX (Fürst et al./2016 2017b). 


Use of the He II 24686 emission line to determine the compact object mass via dynamical means is problematic, 
as{Fabrika et al.|(2015) show that this line in ULXs probably originates in the outer regions of an outflowing wind, as 
in $S433. There remains a question about the width of the He II lines in NGC 7793 P13 — these cannot originate from 
the companion star but also appear to be modulated on the same 64 d timescale as the optical flux 
(Motch et al.[2014); this may point towards the wind precessing as it does in $S433 (see Section 2.3.3 for a discussion 
of super-orbital periods). 


2.5.3. UV 

Models which consider the temperature dependence of the super-critical inflow already im- 
ply that, for high mass accretion rates, the outer photosphere of the wind should emit in the UV around the Eddington 
limit for the compact object. Compton down-scattering of photons from the innermost regions could also produce 
extreme luminosities in the optical/UV band. This may explain the inference that SS433 emits at over 10*° erg/s in 
the UV Waisberg eral [2019 

UV observations of ULXs are limited by the difficulties inherent in observing in much of this band, but a partic- 
ularly notable result is that of who observed NGC 6946 X-1 with HST, finding a spectrally soft 
ULX showing strong variability and QPOs, surrounded by a photoionized and shock-heated optical nebula (MF16 — 
see below). Spectral fitting carried out by these authors implied a UV luminosity in excess of 10°? erg/s, consistent 
with the luminosity needed to produce the observed Hell 24686 emission-line luminosity (Abolmasov et al.|/2008). 

More recently, monitoring of ULXs with Swift has started to provide hints of how the X-rays (via the XRT) and 
UV (via UVOT) may be correlated (or anti-correlated). As shown in {Fuerst et al.|(2021), the ~ 64 d optical period 
observed in NGC 7793 P13 (Motch et al.J2014), is also seen in the X-rays (with a precise value of 65.21+0.15, see 
also{Hu et al.|2017), in the UV (with a period of 63.757017 d), and in the fits to the changing pulse period, yielding an 
orbital period of 64.86+0.19 d. Itis clear from Figure 11 that the UV brightness is not in phase with the X-rays (and 
indeed, the UV modulation is considerably stronger when the X-rays are dim). 


2.6. The environment of ULXs 


ULXs are a product of their host environment, and interact with it through kinetic and radiative feedback. In this 
section we discuss the small and larger scale environments of ULXs and the insights these provide. 


2.6.1. ULX nebulae 

As first reported by {Pakull and Mirioni| (2002), a substantial number of ULXs are observed to be surrounded by a 
bubble of warm gas (although the phenomenology is in reality somewhat varied), visible in the optical via emission 
lines (and in some cases detected in the radio as well). As indicated in Table 1, these ‘ULX nebulae’ are typically 
~10 - 100s pc across (the largest presently known being ~500 pc in diameter) and tend to be considerably larger than 


typical SNRs (e.g. 2014) although there is sometimes disagreement on the actual size (e.g. |Moon et al.[201 1 


24 


800 1000 1200 1400 1600 


Time (d) since MJD 57500 


600 


25 


rIttyittTs TItrtTTTTITTI Jaa el wet = de i] DRA aa ester dee AE eau 
ANE S a 2 erase: J a [2 | A T E: diiic: A EE Ja Jes raS eee J 
E -L L f ag + + 
——— 1 
oa ' 
pe : 
L L L ! L L ENS 
i 
f 
f 
1 i J 
= i 
f 
i 
f 
— ; eo 
i 
Eea O mee NN an E Rr e 
oa 
| TE 
i 
a i 
= ' 
i 
B f a 
i 
' 
=. Ə i 
] ! J 
i 
f 
f 
J i M 
aaan] TE O a AMEN na aaa © MEE R.. 
oi 
i 
1 i L 
i 
f 
f 
i 
i 
a i alr, 
f 
i 
i 
a f 
f 
d ' L 
i 
i 
ee ee Weal | 
J i a i i 
e:i i e| > 
> 1 D> > 1 2 
= 1 1 1 1 
L = L L í L i L i L i 
0 i i i i 
x i i i i 
g= 1 1 1 1 
, € g £ 
g = T g Po T i T 
Ga 1 1 i) 1 
L + L f 4s 4: 4 i 
EE Sieg aero pies ee acct tle ee te des esl ee ee eae 3 
= 7 T 1 an 1 a D 1 TT 1 
— f i l i 
[ ap | i 1 i ale J l 
= 1 t 1 1 
= > ya, ; ~ 4 =~ i i A ~ 
F a BET | oo > v o% p w 
piritiiy potriiitipii ty i bunis tirit ns le ele ee (ieee tet at 
10 a 19 =] 19 r = jmi =| © oo © oo © oo © © 
= nm N © i g 109 wpe ° 
[Seu] LOAN [5/839] over LUX [sn] dV [sn] dv [sn] av [s7] dV 


400 


200 
appear to result in a substantial decrease in the spin-up rate as indicated by the lower panels (see|Fuerst et al.[2021 [for 


the ULX pulsar, NGC 7793 P13, with modulations around a ~65 d period. The drop in X-ray flux in 2019 does not 
details). It is also intriguing that the UV and X-ray modulations appear to be — to some extent — anti-correlated. 


Figure 11 from|Fuerst et al.](2021) showing the long term Swift X-ray (XRT) and UV (UVOT, U band) lightcurve of 


versus|Kaaret et al.[2004). Bearing this point in mind, the dimensions provided in Table 1 do not necessarily account 
for the entire region influenced by the ULX (e.g. there is a reported 800 pc zone of weakly ionised gas surrounding 
NGC 1313 X-1:|Pakull and Mirioni/2002). 

Analogous to the W50 nebula of SS433, we expect gas from the original SNR to be inflated by jets and winds from 
the inner accretion flow (e.g. and be photoionized by the UV—X-ray 
radiation field (which we broadly expect to be anisotropic). Indeed, ULX nebulae show a mixture of photoionized 
(notably He II 46862) and shock ionized lines ([OI] 63004 and [OII] 50072) — some nebulae appearing to be domi- 
nated by one process over the other (although both probably occur in any given source e.g. [Abolmasov et al.[2007a). 
ULX nebulae have conclusively been shown to not be inflated by O star associations (i.e. a superbubble) and are 
unlikely to be due to a hypernova explosion (Pakull and Grisé|2008). 

Applying a model where the ULX nebula is inflated via the mechanical action of outflows (i.e. a shock dominated 
rather than photo-ionised dominated nebula) gives interesting estimates for the age of the system. This in turn provides 
clues as to how the binary may have evolved, and for how long super-critical accretion can be maintained. From 


(1977) the radius of the bubble is (compare Eq. (117)) 
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where R is the radius of the bubble at time f, Liot is the mechanical luminosity of the wind and/or jet, and the ISM 
density (assumed to be uniform) is pọ = pumpn, where u is the mean atomic weight of the gas, mp is the proton mass 
and n is the hydrogen number density. The characteristic age of the bubble is then given by T = 3R/Svexp. 

These equations allowed|Cseh et al.) (2012) to determine the energy content and age of several ULX nebulae, the 
results of which are shown in Table 1. Assuming Ly, = 1 /2MwinaV> im aœ and that the wind provides the entire energy 
reservoir of the bubble (i.e. provides Liot), then equation 18 may alternatively be written in terms of the wind properties 


(e.g. |Pinto et al.!2020a): 


TitwindYorind = 3/5 
a t (19) 


R= 0.76 
Po 


where rwina is the mass loss rate in the wind and vying is the wind velocity. For measured wind velocities in ULXs (e.g. 
[Pinto et al.|2016} |Kosec et al.[2018b), a kinetic luminosity of a few x 10°? erg/s is plausible which gives bubble radii 
~100 pc, consistent with observation (Table 1). We note however that the radius depends on the driving luminosity 
(radiative and mechanical) to the power 1/5, so that even a difference of a factor 100 in luminosity alters the radius by 
only a factor 2.5. 

The He II 46864 recombination line is more easily created via photoionization than by shock-excitation (e.g. 
Berghea et al./2010b), so using this line allows us to treat the nebula as a UV—X-ray bolometer, and infer the irradiating 
luminosity by using photoionization codes such as cLoupy (Ferland et al.| 1998). This has been performed for a number 
of ULXs, notably Ho II X-1 (Pakull and Mirioni|2002|and who used the OTV line in the IR band) 
and NGC 5408 X-1 E ike finding that the required luminosities are a good match to the 
X-ray luminosities observed from the source. 

Such consistency is used as an argument against strong geometrical beaming in these ULXs as the nebula should 
see the isotropic emission (Pakull and Mirioni|2002). However, the ULX wind photosphere emits a luminosity ~ Leaa 
which — depending on the mass transfer rate — can have a temperature ~ 10° K, implying strong near-isotropic soft 
X-rays and EUV (see the discussion of Eqs. (I10|[I1I]below). 

In agreement with this, both Ho II X-1 and NGC 5408 X-1 are spectrally soft and potentially viewed at larger 
inclinations than some of the spectrally harder ULXs. Should this be true then it is not surprising that the irradiating 
luminosity inferred from photoionisation modeling is consistent with the observed luminosity, as we would naturally 
expect the more obscured (at such angles) hard X-ray emission to be the most geometrically beamed (rather than the 
more isotropic soft emission). 

In principle, the number of ULX bubble nebulae can also be used to explore the role collimation/beaming plays 
in the wider population based on the reasoning that, for every ULX observed with a given beaming factor, there 
should be many more which are unbeamed and may not even be detectable as a ULX. As the nebula’s photoionisation 
responds mostly to relatively unbeamed soft X-rays and EUV photons, such a premise would hold only if the nebula 
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. : Inferred age Radio/IR (jet) Lin 
i Dimension(s) (years) j Seine (10°? erg/s) 
NGC 1313 X-1 240 pe 
NGC 1313 X-2 350x500 pc? ~1x10° (b) ~10 
IC 342 X-1 110 pc* ~6x10° Y ~0.7 
Ho II X-1 45 pc? Y 
Ho IX X-1 250 pc®° ~1x10° Y ~10 
M81 X-6 115 x 42 pel 
IC 2574 X-1 
NGC 2403 X-1 300 pc® 
NGC 4559 X-7 
NGC 4631 H7 
NGC 4861 X-1 
NGC 4861 X-2 
NGC 5204 X-1 360 pc 
NGC 5408 X-1 60 pc? ~1x10° (i) Y 0.7 (i) 
NGC 5885 ULX | 200 x 300 pc® 6x10? (j) Y ~20 (j) 
NGC 6946 X-1 20x34 pc* Y 
NGC 7793 S26 | 185 x 350 pc (D) 4x10° (1) Y ~10-100 (l,m) 
SS433 (W50) 50-70 pc” ~2x 10° (b) Y > 0.1 
M83 MQI 13 pe? 1.3-2.0x10* (0) Y ~ 30 (n) 
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Table 1 The properties of ULX nebulae, with non-ULX (i.e. hidden ULX) systems included at the bottom in a 
detached table. Objects and values from|Pakull and Mirioni|(2002)} unless indicated otherwise. Dimension is diameter 


was energised through shocks (which would indeed seem to be the case: |Abolmasov et al.|2007a}|Gurpide et al.!2022) 


or if the accretion rate was high enough such that the photosphere emitted in the UV (with any substantial X-ray 
emission obscured or reprocessed). This then implies we should detect bubble nebulae with X-ray faint (non-ULX) 
counterparts; indeed there is already support for such observations both in other galaxies (e.g. and 
within our own (SS433: see Section 2.8.3). Searches for X-ray quiet bubble nebulae are made somewhat harder as 
ULX nebulae with strong photoionisation signatures (NGC 5408 X-1 and Ho II X-1) appear to be smaller (10s of pc 
rather than 100s of pc: Table 1), precluding detection of such relatively small structures in the absence of high angular 


resolution, narrow filter instruments (Russell et al./2011). However, |Pakull and Grisé| (2008) searched VLT images 
of NGC 1313 for additional bubble nebulae without ULX counterparts and found only one possible candidate. As 


shown by|Wiktorowicz et al.|(2019), should geometric beaming be strongly related to accretion rate (King|2009), then 
binary population synthesis would imply that (in the absence of precession) the difference between the observed and 
underlying ULX population should be a factor of 5-15 unless black holes are the dominant population (in which case 
the observed is very similar to the underlying population). This would imply ~10-30 ULXs should be present in NGC 
1313 (assuming neutron stars dominate the underlying population), around 1/3 of which would have large bubble 
nebulae (Pakull and Grisé[2008). This is not significantly inconsistent with observation and any tension is relaxed if 
the population demographic is not 100% neutron stars or if the beaming is not as extreme a function of accretion rate 
as assumed 

A number of ULXs have radio nebulae (see Table 1) of a similar size to their optical counterparts 
[Cseh et al.[2012). In Ho II X-1 (which shows a clear photo-ionisation dominated optical nebula), the radio 


nebula has been resolved into a discrete lobe-core-lobe structure (Cseh et al./2014). The presence and variability of 
the radio structure in Ho II X-1 implies that repeat ejections also energise the surrounding gas (Cseh et al.|2015) in a 
similar manner to SS433 (W50), the inflated bubble nebula in M83 S2 (Soria et al-[2020), ie) Sa BOT 
and the collimated jet source in NGC 7793 (S26: [Soria et al.[2010). 


2.6.2. ULX host environment 

There is a trend favouring increased numbers of ULXs with lower host metallicity (e.g. and 
higher host star formation rate (see Section 2.1). This makes it important to explore the sub-galactic environment 
where ULXs are located, as this can provide powerful constraints on formation channels. 

Besides the ionised environment (which can be extremely large, e.g. the 800 pc region around NGC 1313 X-1: 


Pakull and Mirioni|2002), it has been widely observed that ULX counterparts tend to be coincident with young stellar 
clusters or OB associations (e.g. |Zezas et al.|2002 2002 2003 2006} |2008} |2012 
[Swartz et al./2009} [Liu et al.|2007} |Abolmasov et al.|2007b). Although ionizing, the stellar environment provides 


insufficient flux to explain the emission line luminosities seen from the bubble nebulae (and the expected kinetic input 


due to supernovae is also not enough to inflate them: 2006) — indeed, in the case of IC 342 X-1, there 
are no O stars within 300 pc of the ULX 

ULXs are also found (albeit rarely) in globular clusters. Relatively few (17 in total) are known at the time of 
writing: three in NGC 1399 (Shih et al./2010} [Irwin et al./2010}|Dage et al.]2019), two in NGC 4649 (Roberts et al. 
2019), five in NGC 4472 (Maccarone et al./2007||2011) and a further seven in Mo Deal 
22020). 

ULXs located in globular clusters are expected to differ markedly from those located in the wider population. 
Given the dense environment of a globular cluster, such ULXs likely have a dynamical origin - 
distinctly unlike those ULXs located in star forming regions (e.g. |Brightman et al./2020b), and the companion stars 
are also considerably older (and therefore less massive, e.g. the ULX may be feeding from a white dwarf: 
potentially in a short period, ultracompact binary [King]2011). In terms of brightness, the typical X-ray 
luminosity is ~10*? erg s~! with the only globular cluster ULX found to exceed 4x10% ergs! being GCU8 in NGC 


1399 (Dage et al. 2019). 


2.7. ‘Hyperluminous’ X-ray sources? 


The term ‘hyperluminous’ has sometimes been used to refer to sources with inferred (isotropic) luminosities in 
excess of 10*! erg s~!. However there is little evidence to support the implicit suggestion that these bright ULXs form 
a physically distinct class. They certainly do not need to contain unusually high mass black holes or IMBH: the ULX 
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pulsar NGC 5907 ULX-1 reaches such luminosities (e.g. and evidently contains a neutron star. 
The well known extreme source SS433 is probably a ULX viewed ‘from the side’ (i.e. from outside the X-ray beam: 
see Section and inferred on evolutionary grounds (see [King et al.]2000a) to be supplied with mass at a rate 
~ 1075 > 10°*Moyr7!. From the beaming formula (55), if viewed from within this beam, SS433 would have an 
inferred isotropic luminosity Lsph ~ 10*'erg s7! if the accretor is a 10 Mo black hole. We discuss the source HLX-1 in 


Section (3.16.2) below. 


2.8. Galactic ULXs 


A number of Galactic sources are known to reach or exceed their Eddington limit for extended periods of time (i.e. 
discounting type I bursts etc). Their proximity allows for more detailed tests of how super-Eddington mass supply 
operates. Here we report those sources which would qualify as ULXs if placed in another Galaxy, with observed (or 
inferred) X-ray luminosities > 10°? ergs™!. 


2.8.1. Low mass X-ray binaries 

Several low mass X-ray binaries (LMXBs) containing black holes are known to exceed their Eddington limit 
during outbursts mediated by the classical disc instability (Lasota|2001), recently modified to account for large discs 
powering super-Eddington outbursts (Hameury and Lasota|2020). GRS 1915+105 reaches rates close to its Eddington 
limit (and has been in outburst for at least 30 years: 1992), GRO J1655-40 may have entered a 
super-Eddington state in its 2005 outburst (Neilsen et al.|2016), V404 Cygni is thought to have exceeded its Eddington 
limit in its recent 2015 outburst (Motta et al./2017| and potentially in the earlier outburst if 
not all of the radiation was visible to us) in a similar fashion to V4641 Sgr in 1999 (Revnivtsev et al.|2002). In addition 
to these confirmed black hole systems, neutron star LMXBs certainly exceed their own Eddington limit and may even 
approach the empirically defined ULX threshold of ~ 10%? erg/s (for instance the bursting pulsar, GRO 1744-28: 
[Sazonov et al.]1997). It is notable that all of these systems have fairly long orbital periods (2.6 - 31 d), which results 
in a large outer disc radius and the creation of a large reservoir of material for accretion (see 
[2020] for details). 


In addition to large X-ray and optical luminosities, and powerful mass-loaded outflows, optically thin, flaring jet 
emission is observed to coincide with outbursts of those sources approaching their Eddington limits 
[1997). Such radio-bright, highly variable events provide another way of detecting extra-galactic systems entering a 
near-Eddington/super-Eddington state (and may favour the detection of black hole over neutron star systems: 


2013}|/van den Eijnden et al.|2018). 


2.8.2. High mass X-ray binaries 

A number of accretion powered pulsars fuelled by high mass donor stars (i.e. high mass X-ray binaries: HMXBs) 
are known to reach x10°" ergs™!, including the X-ray millisecond pulsar AO538-66 (Skinner et al [1982} SMC X-1 
(e.g. |Bonnet-Bidaud and van der Klis} 1981) and SMC X-3 (mentioned already in this review). 

A recent important observation of Galactic HMXBs in the context of ULXs is that of Swift JO243.6+6124 
fet al.|2017). This magnetised neutron star (with an inferred dipole field of > 10!7G:|Doroshenko et al.|2018) accretes 
from the decretion disc of an O9.5Ve star (Reig et al.|2020), reaches a (0.5-10 keV) luminosity of 10°? ergs~!, shows 


pulsations on a timescale of ~ 10s (Jenke and Wilson-Hodge|2017), shows evidence for fast outflows (van den Eijnden 


2019} 2019), and produces radio jets (van den Eijnden et al.|2018). These jets are relatively weak radio 
emitters, with a maximum flux density of 100 uJy at 6 GHz (at an assumed distance of around 7 kpc: |van den Eijnden 


2018), and are therefore unlikely to be detected in other galaxies in the absence of substantial Doppler boosting. 
However, detections of jet events from Galactic HMXBs may become more widespread with the advent of sensitive 


all sky radio monitoring (SKA) and will permit tests of neutron star jet models (e.g. |Parfrey et al.]2016). 


2.8.3. SS433 


Discovered by|Stephenson and Sanduleak|( 1977), the A type supergiant companion star in SS433 (Gies et al 2002) 
is transferring mass to the compact object at a rate ~ 1075 > 1074 Mo yr“! (Shklovsky|1981}/King et al.|2000a 


2006). This would imply ñm ~ 100 and 1000 for a 10 Mo black hole or 1.4 Mg neutron star respectively. This 
strongly super-critical mass transfer rate suggests that $S433 is ‘intrinsically’ a Galactic ULX similar to those with 


29 


luminosities > 10*°ergs~! seen in other galaxies (Begelman et al.|2006). However it does not qualify as a ULX on the 
grounds of observed luminosity — see Section|3.14.1 


One of the most dramatic features of SS433 is the baryon-loaded jet (Fabian and Rees) 1979 1979 
1994} [Marshall et al.|2013) which is ejected at ~0.26c and is twisted into a helical shape by the 162 day 


precession of the surrounding disc and wind. The baryon-loading, low velocity (v/c ~ 0.25), and 162.4 d precession 


period of the jet might result from a collision with the precessing outer disc (Begelman et al.12006) or could simply 


be due to the jet being the highly collimated wind from the innermost regions (which matches simulations, e.g. Jiang 


et al./2014||2019b) tracing the precession. 


We limit ourselves here to a discussion of the similarities and differences between SS433 and other ULXs, and 
direct the reader to for a comprehensive review. 

SS433 can be regarded as a key example of a 
ULX in our own Galaxy, merely inclined such that 
we do not see directly into the wind-cone (Begel- 
[Middleton et al 20216). As indi- 
cated in the previous sections, besides the similari- L owl Ne X OLD (LOW-INT) | 
ties of the broad He II recombination lines n SE 
et al.[2015) to those seen in other ULXs, there are a 
number of sources which look distinctly like $S433, 
viewed at similar edge-on inclinations. Notably the 
inflated nebula in M83 and collimated jet source in 


NGC 7793 (Pakull et al.][2010} [Soria et_al.|/2010 co 
2012). In both of these cases, the jet | NEW (LOW-INT) | 


power appears to be greater than in SS433 for simi- 
lar or slightly higher apparent X-ray luminosities. In- 
deed, very many ULXs must also be oriented in such 
a manner such that we only observe a relative pro- 


portion of the entire ULX population (see|Middleton t ; ; 
and King}2017;|Wiktorowicz et al.|2019). Whilst it is 
L NEW (INT-BRI) 


rather hard to study the innermost regions of S8433 eL ' J 
(although we can study the reflected emission and in- [ ] 
fer the intrinsic flux: [Middleton et al.[2021b), it pro- 
vides a unique opportunity to explore how such sys- 
tems are fed, its feedback into the local ISM, the na- 
ture of the binary (e.g. probes of circumbinary struc- 


ture: |Waisberg et al./2019), and the structure of the 0.5 Boas. as 2 
X-ray obscuring regions (Middleton et al.|202 1b). = 


2.8.4. Future ULX observations . : . . 

It is apparent that the study of ULXs has made Figure 12 F TO (2020b}, showing simula- 
dramatic progress as instrumentation has improved tions which highlight the spectral quality of observations 
(e.g. the discovery of pulsations, winds and jets). In of ULXs (in this case NGC 1313 X-1) which will be made 
the near future, a number of developments will see possible by the launch of Athena with its X-ray integral 
the field evolve in new, exciting ways field unit (X-IFU). The simulations utilise absorption and 

, f emission models derived from existing observations (see 


2.9. X-rays 2020b|for more details). 


Much of our understanding of ULXs derives 
from the X-ray band (even if this only actually grants 
us a restricted view of the full energetics of the sys- 
tem). Several new X-ray satellites are anticipated in the next two decades, including the Astro-H recovery mission, 
XRISM (dedicated to high energy-resolution X-ray imaging: |KRISM Science Team|2020), and ESA’s flagship Athena 
mission (scheduled to launch in the 2030s: [Nandra et al./2013), comprising a wide field imager (WFI) and X-ray 
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integral field unit (X-IFU). Both instruments onboard Athena represent a significant improvement in effective area 
for their respective specialisms: the WFI is around an order of magnitude more sensitive than XMM-Newton’s EPIC 
pn detector (from ~0.5-2 keV), and the X-IFU is between 45 (at 1 keV) and 6 (at 7 keV) times more sensitive even 
than XRISM. This improvement in sensitivity will allow considerably deeper searches for pulsations in known ULXs, 
placing important constraints on the underlying population demographic. The sensitivity will also allow new searches 
for QPOs and reveal new details in the lag spectra 
Middleton et al.|2021b). Finally, the impact on studying ULX outflows and how they vary with spectral state will be 
dramatic, allowing results previously obtained in ~day exposures to be obtained in only a couple of hours (see Figure 
12 which shows a simulation of the atomic lines detected in NGC 1313 X-1 as seen by Athena: [Pinto et al.[2020b). 


2.10. Multi-wavelength 

As we have discussed in the preceding sections, multiwavelength observations of ULXs have been highly reveal- 
ing, allowing the discovery of jets, studies of optical nebulae and likely donor stars, and have provided insights into 
the nature of the accretion flow (irradiation and emission from the outer photosphere). 

The future of optical observing of ULXs — as with other variable astronomical objects of interest — is poised to 
change dramatically with the introduction of all-sky, high angular resolution, high-throughput observing via the Vera 
Rubin LSST, which is due to start operations in 2023. LSST will provide highly sensitive snapshots (e.g. a limiting r- 
band magnitude of 24.7 in a single 30s exposure at SNR = 5) of ~18,000 deg? of the Southern sky (DEC = -65 — +5). 
Combined with a median seeing—limited angular resolution of 0.7”, this will allow deeper probes of possible optical 
counterparts and explore optical variability for signs of irradiation or precession (e.g. 
fet al.[20195). LSST is not equipped with narrow filters (nor any spectroscopic instrumentation), precluding locating 
new bubble nebulae, but the MUSE IFU on the VLT is set to transform our understanding of ULX feedback by 
spatially resolving different regions of ULX bubble nebulae (e.g. |Gurpide et al.|2022). Finally, the introduction of the 
ELT (sometime after 2027), will also provide a step change in signal-to-noise, allowing deep spectroscopic studies of 
ULX counterparts. By identifying lines from the companion star, mass functions will no doubt follow. 

At longer wavelengths, there will be a great deal of interest in the JWST view of ULXs, as this will provide 
opportunities to locate further red super-giant companions (e.g. [Heida et al./2016), study the surrounding nebulae 
(e.g. using the OTV line: at a resolution even higher than can be achieved by MUSE, search 
directly for IR bright ejecta or indirectly search for the presence of jets via IR variability (e.g. |Dudik et al.|2016). 

At the longest wavelengths, SKA will transform our view of the variable radio sky. In relation to ULXs, the 
enormous sensitivity of SKA will allow us to locate LMXB ULXs entering outburst in nearby galaxies (Middleton 
et al.|/2013), search for variability in radio bubble nebulae, and potentially resolve jet structures 
2015). In terms of the impact on studies of Galactic ULXs, SKA will allow the deepest searches for jet ejections 


associated with accreting pulsars (van den Eijnden et al.|2018) and allow models for their production to be tested (e.g. 
Parfrey et al.|2016). 


3. Theory 
3.1. Definitions and notation 
The Eddington luminosity is 
4nGMmyc 33| M E 
Lega = ————— = 1.26 x 10°°| — | ergs, (20) 
o M 


T © 


where m, is the proton mass, or the Thomson scattering cross-section. In this review Lega defined by Eq. isa 
unit of luminosity depending on accretor’s mass only. 
The Eddington accretion rate is defined as 


M 
Mgaa = 1.4 x 10'877} (a) gs! (21) 
M 
=2.2 x 108751 (5) Moyr!, (22) 
© 
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r Spin up rate (ss) . i . Lmax M2 
ULX Ppulse (S) (baseline) Pulse fraction (E range) Pow Py_orb (d) (x10 erg s7!) (Ma) 
ibe ZF i 5-13% (3 - 30 keV) z 
M82 X-2 1.4 1 x 107!! (2.5 years) 8-23% (10 - 30 keV) 25 60 20 >52 
3 E 12% (0.2-2.5 keV) 
_pae 2 10 P 9 
NGC 5907 ULX-1 1.1 8x107! (11 years) 20% (7-30 keV) 5.3 78 200 ? 
NGC 7793 P13" 0.4 4x107" (4 years) 18-22% (0.1-12 keV) 64 64 ~ 10 18-23 (B9) 
i p ~60% (0.2-10 keV) ji 
i, 7 a m 
NGC 300 ULX-1"# 31.6 -6 x107" (3.6 days) ~ 10% (3-20 keV) > 290 ? 5 8-10 (RSB?) 
SMC X-3%1m.n 78 ~-4x10 (110 days) 30-100 % (3-70 keV) ~45 ? ~3 >3.7 (Be?) 
M51 ULX-7°? 2.8 -1x10 (13 years) <5 - 20% (0.1-12 keV) 2 38 ~10 >8 
NGC 1313 X-29 15 1.2x10-7 (98 days) ~5% (0.3-10 keV) <4 ? 20 ? 
NGC 2403 ULX™ 18 3.4x10-% (7) ? 60-100 ? 12 Be? 
Swift J0243.6+6124° 9.86 2.2x10-% (%) ? 28.3 ? >15 Be? 
RXJ0209.6-7427' 93 1.165x10-™ (7) ? >50 ? 1-2 Be 
M51 ULX-8" ? ? ? 8-400 ? 2 400) 
NGC 1313 PULX” 765.6 ? 38% (0.3-7 keV) ? 16 Be? 


Table 2 Key parameters of interest for the ULX pulsars with values taken from those papers indicated by the reference 
attached to the ULX name. The period over which the spin-up/down rate is calculated, is given in parentheses, as is 


the energy range for the measured pulse fraction. a: (2014), b: |B 
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where ņ = 0.170. is the radiative efficiency of accretion, Kes the electron scattering (Thomson) opacity. We will often 
use accretion rate measured in units of Eddington accretion rate and stellar mass in solar masses: 


M M 
and m= —. (23) 


m= — 
Mega Mo 


3.2. Discs with super-Eddington mass supply rates 
At high local accretion rates 7 > 1, the radial structure of a stationary (M = const.) disc around a compact object 


can be divided into three parts (Shakura and Sunyaev}||1973). In the outermost regions the pressure is dominated by 


gas and the electron scattering contribution to the opacity can be neglected. For radii 
Res < 5 X 10°? Rg, (24) 
electron scattering is the main source of opacity, but gas pressure is still dominant. Finally, for radii 
Ror < 1.9 x 10776! “8/71 R,,, (25) 


the pressure is dominated by radiation, while the opacity is mainly due to electron scattering. R, = GM /c? is the 
gravitational radius. 


The standard accretion disc model (Shakura and Sunyaev}|1973 2002) applies only to geometrically 


thin discs, i.e., discs whose height satisfies the condition H « R. When the luminosity approaches the Eddington 
value this approximation breaks down. 
The approximate vertical dynamical equilibrium equation 


— (26) 


where cs = +/P/p, is the isothermal sound speed, P and p are the (total, gas + radiation) pressure and density 
respectively, and vx = ~GM/R, can be written as 


H _ LR) 
T (27) 
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To derive Eq. (27), we have used the expression for the radiative pressure P, = 40-T+/3c and the solution for the local 


vertical radiative transfer equation iii = (3/8)r Th 2016). 


Therefore already at L ~ 0.1 Legg the thin-disc approximation ceases to be valid and for L > Lgaa one has H > R 
so that the accretion flow is more spherical than disc-like. 


The point at which the local flux takes the Eddington value is called the spherization radius (Shakura and Sunyaev 
1973) and can be seen as the boundary separating the thin disc from a more spherical inner accretion flow. 


3.2.1. Spherization radius 
In a stationary disc (M = const), the local radiative flux generated by the disc’s (anomalous) viscosity becomes 
equal to the local Eddington flux at the spherization radius Rspn. 


The viscous flux (Frank et al.|/2002) can be written as 


3 Re n5 Rin 1/2 
Frise = oh Mc h -(5) F (28) 
and the Eddington flux is 
L 
Fea = [ooo (29) 
From Fyisc = FẸEaa, one obtains 
Raph = 15m Rg, (30) 


We have neglected the boundary term in square brackets since Rspn > Rin for rn >> 1 and Rin ~ 6Rg. The spherization 


radius defined in Eq. is 1.1 times larger than the spherization radius in/Shakura and Sunyaev|(1973). 
assume Rsph = (Mc /Leaa) Rs, where Rs = arf] which gives a value of Rsph 1.3 times larger than 
in Eq. (50). 

The outer regions of accretion discs around compact objects can be strongly X-ray self—irradiated. Then the 
condition for Rgpp is 


Feaa = Fin + Frise, (31) 
where we put (see 1999): 
Minc? 
Fin = CO, (32) 


where C is a constant containing all the information about the disc irradiation process and Min is the accretion rate 
onto the accretor. The spherization radius becomes 


in\~ 


1 
Ron =15(1- C7) Re (33) 


For “standard” values of C ~ 107? (Dubus et al.| (2001) disc irradiation can be neglected in the definition of the 


spherization radius for 7 « 200(0.01/C). One should notice, however, that at very high accretion rates it is nor 
clear how self-—irradiation can proceed. Even at sub—Eddington accretion rates when applying irradiation models 
to observed systems, one concludes that scattering in the wind is still not sufficient to produce the observed disc 
heating, even in combination with direct illumination (Tetarenko et al.|[2020). On the other hand, until now, the disc 
considered in this context were rather smaller than those supposed to be present in many ULXs so it this latter case a 
component of order Lgaa emerging from the photosphere of the spherical wind, might irradiate the outer parts of the 
disc. Fortunately, as argued above, this effect is presumably not important at most luminosities of interest. 

High accretion rates implying large optical depth and high radial velocities in the inner regions of the flow can 
result in photons being trapped in the accretion flow. 


6Tn the literature both units of length, Rs and Rg are used. 
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3.2.2. Trapping radius 
The trapping radius is defined as the place where the photon diffusion (escape) time Ht/c equals the viscous infall 
time R/v, . 
Hrv, HK M 
c — c 2mnRÈ 


where È = 2 f pdz is the disc surface density. For R < Rtapp photons trapped in the flow are entrained (advected) 
faster than they can escape, so that if the accreting body is a black hole, most of the accretion energy is accreted by 
the black hole and adds to its mass. In this state there is in principle no upper limit to the accretion rate M but the 
accretion luminosity grows only logarithmically with M (see Eq. rc In other words, very super—-Eddington mass 
supply rates do not imply (very) super—Eddington luminosities. 

Eq. |34/assumes that the accretion flow moves only in the radial direction, and that there are no outflows. 
proposed an alternative to such an advection—dominated solution (which they considered “im- 
probable”) where for R < Rgpn the local flux is everywhere equal to its Eddington value. 


H 
Rerapp = = 20 R mR, (34) 


3.3. Shakura-Sunyaev (1973) second (windy) solution 
To keep Fyis = Fgaa, everywhere for R < Rspn, one must arrange that M ~ R (see Eqs. [28]and[29). Hence 


m(R) = mo (35) 
Rsph 
Using Cr = Fyise = Fgag one can calculate the luminosity of disc described by the windy SS73 solution 
Rsph 
Lihick = a oTi,20RdR x Leda ln m. (36) 
Rin 
The total disc luminosity is then 
Liota =Lihin + Lihick = 

Rsph Roo 
an Í oTi, RAR + f ory Ran] x Lpaa (1 + Inr). (37) 

in sph 


In this SS73 ‘windy’ solution, maintaining the radiation flux at the local Eddington value requires the accretion rate 
to decrease with radius as M ~ R. This requires that matter is expelled from the disc flow as a wind. But the trapping 
radius is close to the spherization radius, so one can also legitimately consider discs in which the “excess” radiation 
is trapped in the inflow and advected on to the accretor, while the accretion rate remains constant. We examine these 
advective solutions in the next two subsections. They have played a very important role in the understanding of 
high-rate accretion onto compact bodies but, as we will see, they do not give suitable models for ULX accretion. 


3.4. Slim discs 


For a black hole accretor, the advected energy disappears behind the event horizon (growing its mass) and we can 
regard advection as an additional cooling mechanism. In contrast, for a neutron star accretor, the advected energy 
must be radiated from the stellar surface or the magnetosphere. Then the accretion luminosity violates the Eddington 
limit and the assumption of steady radial infall. 

The energy conservation equation in a stationary disc is 


Fyisc = F, F Fray, (38) 
where F, is the radiative flux and F,qy is the advective term 2016) 


=xTov, ds M 
ps = oe, 
R dinR aR? (33) 
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Here È is the surface density, T is the midplane temperature, v, the radial velocity, cs = »/P/p the speed of sound, s 
the specific entropy and é, ~ 1 is a slowly varying function related to the entropy gradient. From Eqs. [28] [26] ana [39] 


one obtains i z 
Fay Cs 7 (3) 


F Vise QZR? ~ R 


(40) 


so that for thin (H/R « 1) discs the advective term can be neglected. But at luminosities approaching Legg this term 
become dominant, leading to 
Fisc © Fadv. (41) 


Solutions of this equation represent advection dominated accretion flows (ADAFs), which for high accretion rates are 


called “slim discs” (Abramowicz et al.||1988}|Sadowsk1||2009 2011}|Sadowski et al.)/201 1). 


Solutions representing vertically—integrated slim-disc structure have the form 
m~ Ker! ad (42) 


on the M(£) plane. Kes is the electron-scattering opacity and a the disc viscosity coefficient. 
For these pure advection solutions, one has H/R ~ 1, independent of accretion rate and radius (see Eqs. 


[41] and [Lasota] |2016). But since Fay ~ m? while F, ~ m!’ (Lasota et al.| |2016}, below m < 1, radiative cooling 


overtakes advection and the disc structure switches to the standard radiation-pressure-dominated SS73 disc: 
i ~ kir P ay !. (43) 


Using Eq. the radiative flux of the advection-dominated flow becomes 


L 
= 4 Edd 
F, =0T sy © Tat?” (44) 
so that the luminosity of this part of the accretion flow is 
_ iat 4 Rirans r 
Liit =2 OT 2nRdR x Leda -In x% Leaa ln m. (45) 
Rin in 


Equations and both result from the mechanical equilibrium equation of a radiation-pressure dominated disc 


c c GM 
Er * —- > 
Kes Kes R? 


for H/R ~ 1, and g; ~ (GM/R?)H/R (Paczyński| [1982} [Poutanen et al.|[2007). But although the luminosities are 


similar, in the first case only 7 ~ 1 is accreted on to the compact object while in an advection dominated flow the 
whole of 7 > 1 arrives at the accretor’s surface. 

Since H/R in slim discs is independent of the accretion rate, one cannot get so—called ’thick discs” or ”Polish 
doughnuts” by increasing the rate at which matter is supplied to the compact object. Contrary to assertions found 
in the literature (e.g., [Poutanen et al.| [2007) slim discs are not the same objects as Polish doughnuts. As we have 
seen, they do not significantly increase the true emitted luminosity above Lega, and do not collimate or beam it such 
that the apparent luminosity may be super—Eddington in a restricted solid angle. Accordingly they do not offer an 
explanation for the defining feature of ULXs. However for a black hole ULX they may describe the fate of some of 
the super—Eddington mass supply. 


(46) 


3.5. Polish doughnuts 


Unlike slim discs, for which H/R < 1, Polish doughnuts (Kozlowski et al.| |1978} 1978 
Jaroszynski et al.\}1980}|Paczynsky and Wiita\|1980) have by construction H/R > 1, with long central funnels along 


which most of the radiation is emitted. They were devised in the late 1970s and early 1980s to try to explain the high 
luminosities and jets of QSOs by assuming that these objects are fed at highly super-Eddington rates (see, e.g., 
{1981}. They are often confused with slim discs and incorrectly assumed to be advection dominated accretion flows. 
They have been invoked as solutions to the ULX problem, but we will see that this is not correct. 
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Polish doughnuts are models of stationary and axially symmetric accretion structures around black holes. All 
their properties are obtained from a single function €(R) describing the specific (per unit mass) angular momentum 
distribution at the doughnut’s photosphere. There is no explicit assumption about the doughnut’s interior. One assumes 
only that a.) the photosphere coincides with an equipressure surface. b.) the specific angular momentum at the 
photosphere is assumed to be given by some function €(R) (one often assumes ¢(R) = const.), c.) radiation is emitted 
at the photosphere at the local Eddington flux, that is, the local flux F, is given by 


c 
F, = - Son: (47) 
K 


{Paczyński} |1982),where boldface symbols denote vectors. The specific angular momentum is assumed to be Keple- 
rian at the inner and outer boundaries: £ (Rin) = €x (Rin) and £ (Rout) = €x (Rout), and the inner radius lies between the 
innermost stable circular orbit (SCO) and the innermost bound orbit (IBCO), and Rico < Rin < Risco. Because of 
the strong radial pressure gradients (negligible for thin discs) needed to have H/R > 1, the inner flow edge is pushed 
inside the ISCO. 

The doughnut shape is obtained from 


dH ser) 
— =-|— = f(R, H), 48 
aR G p J(R, H) (48) 


where R,z are the cylindrical coordinates, and the radial component of the effective gravitational acceleration geis a 
function of (R). 
The inner part of the doughnut forms a funnel with an opening angle obeying tana = 1/y ~ ./8p, where p = 


(Rin — Rtpco)/Rs (0 < p < 1) and y = (A/R) max (Wielgus et al.||/2016). 


The luminosity is mainly emitted from the funnel, i.e., 


Lfunnel 


1 
= Linne. 49 
1 — cosa b . p) 


Lapparent ~ 
With b = (1—cos œ) ~ p, the funnel can be very narrow, but this does not lead to an increase of the apparent luminosity. 
Indeed, 


1 Erad v š 
Lapparent X g mnel = = Me x const. Mc’, (50) 


because the accretion radiation efficiency aq ~ p, since by definition the binding energy at Rigco is equal to zero. 

So although Eq. at first sight resembles Eq. (3). in Polish doughnuts one does not have Lapparent >> Leue. This 
means that, contrary to initial appearances, Polish doughnuts are not relevant to BH ULXs (Wielgus et al.|/2016), and 
we already know that they are not applicable to NS ULXs. Together with the similar conclusions we arrived at in 
Subsection for slim discs, this severely narrows the range of potential explanations for ULXs. 

We are left with only one possibility for making an accretion flow appear super—Eddington while simultane- 
ously ensuring that the mass supply rate very close to the accretor is not significantly super—-Eddington. In the next 
Subsection we ask if the disc winds required by the Shakura—Sunyaev ‘windy’ solution can collimate the modestly 
super—Eddington intrinsic luminosity Lgga(1 + Inm) (cf given by a super—Eddington mass supply. We note that 
Eq. is an approximation since the spherization radius is not a uniquely defined quantity, which we should not see 
as a rigid limit defining the boundary of the wind. Eq. is valid in both the SS73 and advection—dominated cases, 
but when these effects are simultaneously present one gets a slightly different formula (see e.g.,/Poutanen et al.|2007). 
In the following we therefore use L = Lgaa(l + Inm) as a “universal” reference value for unbeamed luminosity. 


3.6. Beaming by accretion disc winds 


We have seen that in the SS73 ‘windy’ solution, the accretion flow carries only the local Eddington mass rate at 
each accretion disc radius, while the excess is blown away in winds at each radius near the disc centre. These winds 
must be quasispherical, but leave open funnels around the central disc axis (since they cannot achieve zero angular 
momentum). The funnels offer a suitable way of collimating the accretion luminosity, as all other photon escape 
routes have high optical depth. We will see in Section that numerical solutions of such flows are still some way 
from providing a picture which is easy to apply to observed ULXs. 
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Fortunately, soft X-ray observations of ULXs provide a remarkable insight into the beaming by disc winds. The 
soft components of ULXs can be fitted with blackbody spectra, giving the temperature T and apparent (assumed 


isotropic) luminosity Lyon. Kajava and Poutanen) (2009) show that the brightest soft X-ray components (above 3 x 


10% erg s7!) in 9 ULXs cluster around the relation 
Leon = 7 X10" ToT OCS (51) 


where To. key is the blackbody temperature expressed in units of 0.1 keV/k, where k is Boltzmann’s constant. How- 


ever, not all ULXs show such a relationship (cf 2020}|Gurpide et al./2021) and — even in the absence of 


additional spectral components in the presence of a neutron star accretor — the combination of anisotropy/beaming, 
precession and the lack of a suitable model for extracting L,of. unambiguously (see, e.g. [Robba et al./2021), must 
result in complexity and a range of observed slopes. 

At first sight seems paradoxical, as one would expect a blackbody source of fixed size to vary as L œ T4 
instead’ But the relation would instead imply that the radius of the blackbody source varies along with the 
luminosity and temperature. This is of course what we should expect for emission from the central part of accretion 
discs fed at super—Eddington rates, progressively ejecting the excess accreting matter from inside the spherization 
radius, as the latter varies œ 71 (see Section. 1p. We can see this explicitly by parametrizing the true emitted luminosity 
(pre—collimation) as 

L = 4npR’°oT* (52) 


where R = rR, œ M with r = 15m (see Eq. Bo}, and o is the Stefan-Boltzmann constant. Setting L = [Lea « M, 
where / ~ 1, and p allows for geometrical projection, we can write 


2 
rp 
L œ R°T*p « M’T*r’ px PT (53) 
where we use R œ rR, at the first step and M œ L/L at the second. The final form implies L œ T~*. An observer 
assuming that the flux is isotropic with the observed value, rather than collimated, now infers a total blackbody 


luminosity Lsph = b“'L, so that, inserting the constants, 


4n 1 P 28x10“ / P = 
Lh = pel ag = dk z | ees E (54) 
ox T* pbr Tõikey \POr 


The fact that observation gives L œ T~* (cf Eq. 51) means that the factor /?/pbr? in this equation must be a constant. 
Since r = (15/2) (Eq. 30), consistency with Eq. (51) requires 
73 
b= —, (55) 
m 
where we have assumed that the slowly—varying factors l ~ p ~ 1. (Tight beaming means that we observe the source 
along the disc axis, so we view the central disc plane orthogonally, making p ~ 1.) 
The physical origin of the simple form is straightforward. The beaming factor is the total solid angle of the 
two open cones around the disc axis, i.e 
b = (1 — cos 6) (56) 


where 8 is the opening half—angle. Near the spherization radius, where the wind outflow is maximal (see Section 3.3), 
the flow geometry scales with 7. But all flows, regardless of m, reduce locally to the same Eddington inflow near the 
accretor, so we expect @ ~ constant x 7n!. For sufficiently large ri we have 8 < 1 and so b ~ 67/2 œ mn. 

We see that the simple form formally applies only for m = 8, since, at such rates, b < 1. Evidently it is 
legitimate to adopt a form of b interpolating smoothly between b = 1 for m = 0, and the form for mn > V73 = 8. 


For example,|Hameury and Lasota|(2020) use 
73 


ie 7 
T3 +m? 67 


TIndeed one does find L œ T+ for sources which are permanently below 3 x 10° ergs™! (Kajava and Poutanen{2009}. It is likely that these 


sources are black holes with masses large enough that their luminosities are sub—Eddington. 
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Importantly, although we have derived this relation by considering only the soft X-ray emission of ULXs, the fact 
that it is purely geometrical and involves only electron scattering means that it must apply to all forms of radiation 
from ULXs. 

The true emission from the central accretion disc is L ~ Lega(1 + 1n rn), so the apparent luminosity is 


Leda . M 
Lsph = oy. Mi? ~~ M’ (58) 


where M is the mass supply rate, since m = M/Mgaa œ M/M?. So for a given mass supply rate, accretors of lower 
mass have higher apparent luminosities, as they are more super—Eddington, and the consequent increase in beaming 
outweighs the lower Eddington luminosity. This means that neutron stars are favoured as ULXs over black hole 
accretors in binaries where the mass transfer rate is insensitive to the accretor mass. This condition holds for mass 
transfer driven by the evolutionary expansion of the donor star, (see Sections[3.1]and[3. 13), and also for cases where 
the increased mass transfer is a transient effect caused by a disc instability, e.g. in Be—X-ray binaries, and in soft 
X-ray transients, where the disc is subject to the thermal—viscous instability (cf|Hameury and Lasota|2020). 

The beaming formula gives a simple explanation of the observed soft X-ray correlation (51). We can easily 
see that other explanations are problematic. 

If we assume there is no beaming, i.e. b = 1, then (54) requires 


l 


The true accretion luminosity is L = / Lega, and so (54) would require an accreting mass 


500 Mo G) 


4 
To. keV l 


(60) 


Maceretor = 


for each object. 

Without beaming, the typical size, r, of the soft X-ray emitting region cannot self—consistently be much larger 
than O(1) for an accretion disc spectrum. Then since p ~ 1, requires l ~ 1/62, we find Maccretor Z 3 X 104 Mo 
from (60). i.e. a significantly massive IMBH. 

A magnetar model would require a neutron star mass, so requires / > 350. But then from there must 
be an implausibly large soft X-ray region r ~ 1.6 x 104, which is ~ 500 neutron star radif] In addition, the implied 
luminosity is L ~ 5x 10*° erg s71, beyond the domain of validity for magnetar models for ULXs (Canuto et al., 1971). 

In the context of the model described above, the observed correlation (51) provides strong evidence in favour 
of accretion—disc wind beaming (and indeed, yields sensible results: Re Cen We also note that the 
innermost (hottest) regions will naturally be the most collimated and therefore beamed, whilst one would expect radi- 
ation produced from around rsp, to be more isotropic. The final luminosity we observe may therefore be a somewhat 
complicated function and will naturally be inclination dependent. However, regardless of the final functional form of 
the beaming, it is inevitable that where the wind is optically thick and subtends a large scale-height (as seen in every 
MHD simulation of super-critical accretion), collimation and beaming must result. 


3.7. Models of pulsing ULXs (PULXs) 


The discovery (beginning in 2014) of a small group of ULXs showing coherent pulsing has had a significant effect 
on the study of ULXs in general. First, it subjected existing models of ULXs to a stringent new test, as the mass of 
the accretor must be close to 1.4 Mọ. As remarked above, it removed the main motivation for the IMBH model, and 
we will see in Section[3.7. I|below, that it fitted naturally and self—consistently into the the disc-wind beaming picture, 
even though this was not developed with PULXs in mind. 

But the removal of one candidate model (IMBH) did not simplify the discussion, as the discovery of PULXs 
naturally stimulated suggestions that magnetic effects on the accretion flow might be significant. 


8We will see in Section 3.14) that a large photosphere is possible in ULXs producing strong outflows, i.e. if a strong accretion—disc wind 
produces beaming. 
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3.7.1. The KLK17 model of PULXs 

The discovery of PULXs has shown that ULX apparent luminosities can be highly super-Eddington. The lumi- 
nosities of the classic (non Be-X) PULXs range from ~ 20 to ~ 385 times the Eddington value for a neutron star 
(see Table p). For normal (pulsar-like) magnetic fields ~ 10!! — 10!3G, such values exclude the possibility that they 
correspond to the actual luminosities since this would require implausibly high accretion rates (see Eqs. [36]and [45). 
We will show in Sect 3.9] that the hypothesis that PULXs contain magnetars is extremely unlikely. Then the only 
physically reasonable option left is that the luminosity of magnetic ULXs is beamed according to Eq.(58). Given this, 
(2016) proposed a beamed-emission model for the first detected PULX, M82 ULX-2. This model 
was successfully used by 2017) (hereafter KLK17) to explain the properties of the three first-discovered 
PULXs, and later King and Lasota] (2019) applied this model to the growing number of observed PULXs, including 
the Be-PULXs. (2020) explained how one obtains significant pulsed X-ray luminosity fraction in 
beamed radiation (see Section 3.7.5). It is assumed throughout that for neutron star magnetic fields of the standard 
X-ray binary strength (i.e. B < 10! G) the surface emission is locally isotropic. 

noted that PULXs are sharply distinct from other X-ray pulsars, not simply in their 
luminosities (by definition L > 10°%erg s~!) but also in their very large spinup rates (v > 107!°s~? — up to two orders of 
magnitude larger than normal XRPs). Between superoutbursts, transient Be-PULXs are normal spinning-down X-ray 
pulsars, but when their luminosity reaches ~ 10°° erg s~! they become rapidly spinning-up sources with ý > 107!°s~? 
(see Table p}, demonstrating that such high spin-up rates are a generic property of PULXĖ] 
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Figure 13 The L39 — ¥_19 diagram for XRPs and PULXs. Red dots: the ten PULXs with known spin-up rates (values 
from KLK17). Blue diamonds: selected (for comparison) X-ray pulsars (see Table|Appendix Af- in the Appendix) 


The luminosities and spin-up rates of PULXs are tightly correlated (Fig. [13). This correlation strongly implies that 
the dominating torque in the system is provided by accretion, as assumed in KLK17 and|Vasilopoulos et al.|(2018), 


since 
__ J(Ru) _ M(GMRy)'”” 
ii nl 2nI 
where Ry œ M~?’ (e.g. Frank et al., 2002) is the magnetospheric radius and 7 the neutron star’s moment of inertia. 
If we define PULXs by the twin properties of luminosities larger than 10°? ergs! and spin-up rates larger than 
~ 107!°s~?, these two quantities should form the basis for models of these sources. At the very least, any cogent 


oc MS!7 (61) 


°One might insist that it is a coincidence that Be-X systems show the same spin-up as other PULXs when they reach the same luminosities, but 


here one rather expects that “coincidences mean you're on the right path” A 
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model must be able to reproduce these values. KLK17 chose the first option, using the observed quantities Ly, v as 
input for their PULX model. Although there are no direct mass measurements for neutron stars in PULXs, the likely 
mass range is limited. KLK17 make no assumptions about the magnetic field strengths, and these are outputs of the 
model. The resulting values lie in the typical range for standard pulsing X-ray binaries (see Table B}. and compare 
favourably with the dipole magnetic field strengths measured from CRSFs (one of them for a non-pulsing neutron star 
ULX - see also section 2.2.3). 

For the mass inflow, KLK17 use the super-critical (‘windy’) solution of[Shakura and Sunyaev|{1973) described in 
Section 3.3] In this picture, the super—Eddington part of the mass supply below Rspn is eventually expelled as a wind. 
This simultaneously provides the beaming, making PULXs appear super—Eddington and explains the strong outflows 
observed from PULXs. The neutron star gains mass only at its effective Eddington rate (Eq. [36). 

The spherization radius defined in Eq. can be written as 


Rspn = 2.3 X 10° ing; cm, (62) 


where Mo = tito Mgaa is the accretion rate at Rsph, assumed equal to the mass transfer rate. For R > Rgpp the disc is 
assumed to be a standard SS73 accretion disc. From Eq. we have 


F 5 R 
M(R) = hoMeaip—: (63) 
sph 


for R < Rsph- 


KLK17 assume that the [Shakura and Sunyaev| (1973) ‘windy’ model describes the accretion flow between Rsph 
and the magnetospheric radius Ry defined by the equation (Frank et al.|/2002) 


Ry = 1.2 x 108g w” m KN cm, (64) 
where q ~ 1 is a factor taking into account the geometry of the accretion flow at the magnetosphere. Then from Eq. 


(63) 


M 


: _ R 
M(Rm) = Mor 
sph 


(65) 


Applying the model to PULXs with measured spin-up gives Ry < Rspp in all cases. This implies that in these objects 
the accretion flow inside the magnetosphere is highly super—Eddington. Since it is hypersonic but forced to follow 
the magnetic fieldlines it is highly dissipative, and one expects it to generate an outflow similar to that of|Shakura and 


(1973), limiting the local luminosity to its Eddington value. 
The total luminosity is then given by Eq. 


L = Lega [1 + Inmo]. (66) 


The luminosity from both parts of the super—-Eddington outflow is assumed to be beamed. Outside the magneto- 
sphere the beaming factor is taken to be 


pa (67) 
m 
as in Eq. (55). 
For accretion rates such that radiation is geometrically beamed as described by Eqs. and (67), one can deduce 
rno from the observed X-ray luminosity L = Ly by combining these two equations into 


Ly = 2.0 x 10°°rng [1 + In ño] m erg s™'. (68) 


Having ring one obtains R,pn from Eq. (62). 
The second observed quantity, the spinup, follows from Eqs. [61]and|64]as 


PSS 10 ghee (69) 
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Table 3 KLK17 model: derived properties of neutron star ULXs. 


Name tio [|b ||Bqg m P 73? R3 [G]|[Repnm; e| Rum PR [em]||Peqq 6m” [s]]]teq [yr]! 
M82 ULX2 46 110.03 4.0 x 10% 1.1 x 108 1.0 x 107 0.02 15600 
NGC 7793 P13 25 110.12 L1x 10" 5.8 x 10’ 1.6 x 107 0.09 1386 
NGC5907 ULX1 95 |/0.01 9.4 x 10” 2.2 x 108 1.1 x 108 1.86 0 
NGC300 ULX1 24 10.13 5.3 x 10" 5.5 x 107 3.2 x 107 0.19 297 
M51 ULX74 28 [10.09 1.9x 10" 6.4 x 107 2.0 x 107 0.08 1337 
M51 ULX7" 28 110.09 6.9 x 10° 6.4 x 107 4.6 x 10° 0.01 ~ 10° 
NGC 1313 X-2 46 |/0.03 5.3 x 10 1.1 x 108 1.8 x 10° 0.03 8641 
SMC X-3®° 18 110.23 2.3 x 10 4.1 x 107 7.1 x 106 0.006 76621 
NGC 2403 ULX® 13 [10.43 2.5 x 10" 3.0 x 107 2.3 x 107 0.16 578 
Swift JO243.6+6124®*||| 14 110.37 1.3 x 10! 3.2 x 10” 1.7x 10" 0.07 2047 
RXJ0209.6-7427% 17 10.25 5.3 x 10 3.2 x 10’ 1.8 x 10° 0.03 8665 
NGC 1313 ULX8e* 15 110.32 7 3.5 x 10" ? ? 7 
M51 ULX8 17 [10.25 ~3x10!* 3.2 x 10" 2.7 x 107° non-pulsing 


— Systems with ®° superscript have Bestar companions. 
' _ calculated using the value of Pegq77/°m 1/3 from the previous column. 


— consistent with observations 1.)(2018c). 
*_ from observations 8||Middleton et al. [202 1b}. 


— for a ~ 10!? G magnetic field. 


a _ for ý = 2.8 x 10719; ? — for ý = 3.1 x 107!! (Vasilopoulos et al.||2019). 


€ — unknown v. 


This gives the values of the magnetic moment u = BR? (where B is the field and R the neutron-star radius) which in 
turn allows one to calculate the values of Ry and m(R m) from Eqs and (65): 


= 0.04 q7 mi P RP Gem? (70) 
Rm = 1.0 x 103? mi P EL cm (71) 
rn(Ru) = 9.9 om OEE, (72) 


where 107!°y_19 = v. 


The results are presented in Table [3| with the values of the equilibrium period Peg = 0.23q”/ êm Me gine s, and the 
lower limit on the time to reach equilibrium at the present spin-up rate 


fq = Abe > 5) (73) 


(The corotation radius Reo = (GMP? jare)’ is always larger than the magnetospheric and spherization radii. Thus 
only possibly NGC 5907 ULX1 is close to its equilibrium spin, and none of the observed PULXs is a propeller.) 

For NGC300 ULX-1, the KLK17 model predicts a neutron-star magnetic field of B > 5 x 10!?G (since q < 1 
and mı > 1, the entries provided in the 4th column of Table [3]are a factor of few lower than the predicted magnetic 
field values). A magnetic field of very similar strength is deduced by from the CRSF inferred 
to be present in the X-ray spectrum of this PULX. According to the KLK17 model, all seven PULXs have magnetic 
fields in the range 10!! — 10'°G, i.e. always in the standard pulsing X-ray binary range and below the value defining 
magnetars. 

The predicted values of the beaming factor b are in the range ~ 0.03 — 0.5, except for the very luminous source 
in NGC 5709, for which b ~ 0.01. For NGC 300 ULX-1, the model gives b = 0.13: use the 
KLK17 model and obtain b ~ 0.25, because they deduce rno from the average rather than the maximum luminosity. 
The beaming factor b ~ 0.12 for P13 in NGC7793 is consistent with observations of the X-ray irradiation of the 
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neutron-star companion (Motch]|2018). According to the KLK17 model, the luminosity of most of the known PULXs 
is only mildly beamed. This is worth stressing since critics of the KLK17 model often claim that there are problems 
with ‘strong beaming’. 

We see that for all seven PULXs with known spin-up rates ý, one has Rm < Rspn, which is probably the condition 
for observing pulses at all if the mass transfer is super—Eddington. The small difference between R,,, and Ry means 
that the flow is strongly super—Eddington on reaching Rm. Most of this cannot land on the neutron star (let alone its 
polecaps only) and so must be ejected. 

It is worth stressing that nothing in the assumptions of the model guarantees its self—consistency, i.e. that Rsph > 
Ry. But in every case this is satisfied by the deduced values. Further, exactly the same equations describe both the 
Be-star magnetic accretors and those PULXs which do not have Be-star companions: the first group cannot of course 
be modelled at all by assuming very strong dipole magnetic fields. 

The spin-up time-scales teq are probably much shorter than the lifetimes of the individual PULXs (see Table 
for NGC 1313 X-2) so it seems very unlikely that we observe these systems during their only approach to spin 
equilibrium. Instead, although most of them are probably rather far from their Peq they have nevertheless the time 
to alternate spin-up and spin-down phases. Evidently, we can only see these systems during spin-up phases (so that 
ý has its maximum value) because centrifugal repulsion during spin-down presumably reduces the accretion rate and 
so the luminosity. In addition the magnetic fields of the PULXs appear to be significantly lower than is usual for a 
new-born neutron star (10!7 — 10'°G) which is consistent with hypothesis that accretion of even a relatively small 
mass severely reduces the surface fields of neutron stars — this is central to the concept of pulsar recycling, which 
is implicated in the production of millisecond pulsars (see Bhattacharya and van den Heuvel] and references 
therein). Things are even more complicated by the fact that the vast majority of ULXs do not pulse, despite containing 
(magnetized) neutron stars, which suggests that alignment of spin and central disc axes is rapidly suppressing (perhaps 
temporarily?) the pulsations, which are observed to be transient anyway. As expected from the KLK17 model, Be-star 
systems are prominent among the PULXs, since accretion is relatively weak and transient. 

Parametrizing Ry = fRspn, with f ~ 0.3 — 0.9, then from Eq. (65), m(Ry) = fmo (KLK17) and from Eqs and 
(72) we find 

VET Pg em T, (74) 


in agreement with the observed spin-up values of PULXs (cf Fig|13). 

It has been claimed (e.g. that the KLK17 picture suffers from fundamental weaknesses 
since ‘strong beaming’ would be incompatible with the high observed pulsed fraction of the ULX radiation. As we 
noted (two paragraphs above), the beaming is not ‘strong’. Further, in Sect. (3.7.5) using the results of 
(2020), we will show that the KLK17 model is fully compatible with the observed X-ray pulse fraction (which 
varies in time and frequency), and also explains the very small number of observed PULX. 


3.7.2. The Guirpide et al. PULX models 

suggest an alternative scenario in which the hardest ULXs are powered by strongly mag- 
netised neutron stars, so that the high-energy emission is dominated by the hard direct emission from the accretion 
column. It is assumed that one can explain softer sources as weakly magnetised neutron stars or black holes, in which 
the presence of outflows naturally explains their softer spectra through Compton down-scattering and their spectral 
transitions. Outflows are also invoked to explain the dilution of the pulsed emission in sources containing neutron 
stars. According to these authors NGC 7793 P13 and NGC 1313 X-2 would be strongly magnetized with Ry > Rsph 
in contradiction with the prediction of the KLK17 model. However, they find that NGC 5907 ULX-1 appears harder 
when dimmer. This is difficult to reconcile with their scenario and they have to conclude that for this source Ry ~ Rsph, 
in agreement with the KLK17 model prediction. However, (2021) do not take into account the possi- 
bility that the flow geometry is as described in Sect. In this case, under appropriate conditions one sees the 
“magnetic column” directly when Ry > Rgpn. In any case, one would expect that the flux—hardness relation would 
also depend on the intrinsic properties of the emitter and not on the source geometry only. The model geometry of 
is used by|Pintore et al.|(2021) with additional ingredients to interpret the spectral variation 
during the flares of NGC 4559 X7. 

Although M51 ULX-8 is not a PULX, we know that it is a magnetic system. Assuming that its X-ray radiation is 
beamed, one can obtain its rp and so both Ry and Rsph. As in PULXs this gives Ru S Rsph (see Table B}. 
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(see also Sect. speculate that M51 ULX-8 might be a PULX, and suggest that pulsations may be 
seen in future observations of this system. point out that rather long XMM-Newton observa- 
tions would be needed to detect ~ 1s pulsations if the pulsed fraction is < 45%, as in most PULXs. (Interestingly the 
exception is NGC 300 ULX1, in which inferred the presence of a cyclotron feature). 

Although simple, the KLK17 model correctly reproduces the main observed properties of PULXs and has passed 
this test with each newly discovered PULX, in particular with the Be-PULX systems (see Table B). In the KLK17 
picture, PULXs are “normal” (non-magnetar) X-ray pulsars caught in a high accretion-rate episode of their binary 
evolution (non—Be systems: see Sections 3.11, 3.12) or orbital history (Be systems: see Section 3.15). 


37:3: The|Erkut et al.|(2020) PULX models 
(hereafter ETEA20) base their models of PULX on the well-known approach to the disc- 
magnetosphere interaction best known from the paper of|Ghosh and Lamb|(1979). They assume that the luminosity in 
the magnetosphere is at most critical, using, as the relation between the observed X-ray luminosity and the accretion 
rate: 
Ly = TMe, (75) 


where M, is the accretion rate at the neutron-star surface and b is a beaming factor, and assumes M, = Min, where 


My is the accretion rate at the innermost disc radius Rin. (This is separated from the magnetosphere by a boundary 


region of width AR = ôR - b.) 
In other words 
B\3 
1+311 Fa 
B. 
where B, = mc [he = 4.4 x 10°G (see Eq. pa}. 
The beaming considered by ETEA20 is assumed magnetic (i.e. not due to an outflow in origin, with beaming 
factor 


Ix < Lk, = Lg, (76) 


(17) 


where A, is the polar cap area, A = 47R7, is the total surface area of the neutron star, and y < 1 “is a normalization 
constant corresponding to the maximum fractional area of the polar cap” (Erkut et al.|/2020). Since they are physically 
different, and act in different parts of the accretion flow in PULXs, direct comparison of the values of ETEA20’s bm 
with those of KLK17’s outflow-—generated b is not of much interest. 

The accretion flow between Rin and Rspn (defined as in Eq. (30), with a factor with 27/2, instead of 15) is allowed 
to be super—Eddington and is assumed to be described by the ‘windy’ solution (Eq. [63}. 


The inner radius of the disk is determined by the balance between magnetic and material stresses 


d IRZ = 2 p+ 

~ (MR Q) = —R°BSB., (78) 
where Q is the angular velocity of the innermost disc matter within the boundary region, B, is the poloidal magnetic 
field, B; = yB; is the toroidal magnetic field at the surface of the disc, and yg is the azimuthal pitch that can be 
expressed as yg ~ Wy — 1, where w, is the so-called “fastness parameter” 


OQ, Ra \ I 
W, = ($) i (79) 


Qk (Rin) ~ Reo 


At spin equilibrium, the fastness parameter is assumed to be equal to a (model-dependent) critical value wy, = we. 
ETEA20 use we = 0.75. 
Integrating Eq. (78), using Eq.(75) and defining poloidal field B, ~ —u/R? one obtains an expression for the inner 


disc radius i 
PE E a í 


80 
bMRxLx om 
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Since one can write for a strong enough magnetic field 


Bo Sh cos? a, (81) 


A 4Rin 
where a is the angle between the rotation and magnetic axes (Frank et al.|/2002| the star’s spin and the disc rotation 
axes are assumed to be aligned), one gets b = (cos? a/4y) R,/Rin. 


Using u = B-R3 one can then write 


2 | cos’ aLx i (82) 
32y7/2 N\GM,B?R?P6) | 
for the beaming factor and express the fastness parameter as 
yVGM,BRt5\" 27 m 
wx = > 
cos? aLy PVGM, 


where P = 27/Q,. In addition to observed quantities, Ly (ETEA20 use the observed flux Fy, but since, with one 
exception, PULXs are extragalactic, this precaution does not seem to be necessary), and P, Eqs/ [82] and [83] contain 
unobservable (at least in PULXs) quantities œ and y. They can be eliminated by the use of the torque equation that 


ETEA20 write under the form 
27l v = w Pno VGMyReoMs, (84) 


where no is a constant of order unity corresponding to a dimensionless torque. Given the neutron-star mass (ETEA20 
try M, = 2 Mo for M51 ULX-7), moment of inertia and radius, one then obtains two equations for the four quantities 
defining the physics of the system: B, by, wx and 6 as a function of the observables Ly, P and P (or v): 


2GM,1|P 
B= = x! | (85) 
R} nngoô 
and 5 
ila 
= 16n*GM, |_———__—_—_] . 86 
Ox vig * (| ( ) 


These two equations allow us only for delimiting the range of magnetic-field strength and fastness-parameter values 
compatible with observation. The boundary region parameter is assumed to lie in the range 0.01 < 6 < 0.3. The 
determined range of physical parameters can be made narrow by assuming that Ly = Le, or Ly < Le and the fastness 
parameter can be set equal to the equilibrium value wy = we = 0.75. 

In their paper, ETEA20 consider several alternative scenarios: (i) the PULXs are away from spin equilibrium; an 
efficient standard spin-up torque is used to account for the observed spin-up rates, (ii) the PULXs are so close to spin 
equilibrium that the fastness parameter is given by its critical value, (iii) the X-ray luminosities of the PULXs can 
be well represented by the maximum critical luminosity, (iv) the conditions described in (i) and (iii) both apply, (v) 
the conditions described in (ii) and (iii) both apply, and (vi) the X-ray luminosities of the PULXs are subcritical and 
either the condition described in (i) or the condition described in (ii) applies. They found that the narrowest ranges 
for B, bm, and w, are found when they use the critical luminosity condition along with either the observed spin-up 
rates or the spin-equilibrium condition. Scenario (iv), based on the observed spin-up rates at critical luminosity, is the 
only one that works well for all the PULXs, implying B in the 10!! — 10"? G range. The results for scenario (iv) are 
presented in Table [4|(Table 5 in [Erkut et al.|/2020). The ETEA20 models do not allow for a prediction of the pulsed 
fraction of the emitted radiation, nor the pulse shapes, and the authors refer to the “accretion curtain” of 
as a possible source of the pulsed luminosity. 

The conclusions of the [Erkut et al.) (2020) paper: PULXs have sub-critical magnetic fields and their luminosity 
is subject to medium beaming are the same as those obtained by [King and Lasota] (2016); [King et al.| (2017); [King] 
but they differ in their physical motivations. First, the beaming mechanisms are different in 


the two approaches to the problem. Second, while ETEA20 require the accretion flow to be always at most critical 
and assume that the accretion rate is constant inside the magnetosphere, KLK17 allow the accretion flow to be super- 
Eddington and losing mass also when following the magnetic field-lines. It is perhaps worth noting that these two 
very different models both exclude magnetar-strength fields. 
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Table 4 ETEA20 models of PULXs [scenario (iv)] 


Source Lx (erg s71) v Ss bm Ws B(G) 

M82 X-2 6.5 x 10° J1.1.x 10-™] 0.16-0.25 | 0.060-0.22 ]2.0-2.9 x10” 

ULX NGC 7793 P13 || 9.4 x 10°? | 1.2 x 107!! | 0.021-0.027 | 0.065-0.14 | 1.2-3.2 x10!! 

ULX NGC 5907 1.5x 10" | 3.9x 107° | 0.19-0.30 | 0.16-0.61 |2.7-3.8 x10" 

NGC 300 ULX1 4.3 x 10? |5.6 x 107! ~ 1 0.015 6.8 x10!” 

M51 ULX-7 7.0 x 10°? |3.1x 107!! | 0.052-0.0 |0.019-0.063 | 6.2-10 x10!” 

Swift JO243.6+6124 || 1.7 x 10° |2.4 x 107!° ~] 0.041 3.5x10!2 
Note. In this scenario the magnetic field, beaming fraction, and fastness parameter for Swift JO243.6+6124 are only 
marginally obtained at by ~ 1 for M, = 1.2Mo and R, = 10 km, and no common solution is found for NGC 1313 X-2 


(even a mass as low as 0.9 Mo does not help). For all other sources M, = 1.4 Mo and R, = 10 km are used. 


Rsp N 


Figure 14 Schematic picture of a PULX according to|Mushtukov et al.|(2017||2019). Rm and Rsp are respectively the 


magnetospheric and spherization radius. 4 is the angle between the disc plane and the magnetic axis. Accreting matter 
not lost in the wind is accreted onto the central object forming an envelope which is optically thick at accretion rates 
typical of PULXs. In the KLK17 model, the configuration is similar, but Rm < Rsp and the wind does not end up at 
the magnetospheric radius. (see Fig. [15). Note that the Figure of Mushtukov et al., (2017, 2019) shown here assumes 
that the neutron star spin is aligned with the disc axis, although this is not explicitly stated. We argue in subsection 
that this assumption is in general not justified. In addition the configuration shown appears to be axisymmetric, 
i.e. has the magnetic axis parallel to the spin (and disc) axes. This would of course preclude pulsing, so there must in 
reality be some deviation from axisymmetry. 


3.7.4. The origin of “sinusoidal” pulses in PULXs 

X-ray pulse light curves for PULXs are generally described as being ‘sinusoidal’, i.e. without obvious eclipses, 
and continuously modulated (see Sect. |2.3.2). In other words, the X-ray light curves of PULXs are never flat and 
never zero. 


Mushtukov et al.|(2017||2019) argued that this was understandable if, perhaps because of the large optical depth, 


the photospheric surface emitting the X-rays occupied a large fraction of the magnetosphere. Then this region would 


neither be permanently in view, nor periodically occulted. Accordingly, |Mushtukov et al.}(2017||2019) proposed that 


the X-rays are emitted by an optically thick envelope defined by the accretion flow over the neutron—star magneto- 
sphere. This is consistent with the accretion along fieldlines being a bending hypersonic flow, which must therefore 


shock as suggested by|King and Lasota|(2019) (see also Sect|3.7.1). The spin modulation occurs because the magnetic 


axis is inclined to the spin axis. 


Mushtukov et al.| 2017} |2019) do not attempt to reproduce the observed PULX lightcurves and do not explain 


why magnetized systems such as ULX-8 in M51 do not pulse. The simplest explanation would be that the spin and 


45 


magnetic axes are aligned, but then one would have to explain why the few pulsing systems are different in this respect 
from the probable majority of NULXs. 


3.7.5. Beamed pulsed X-ray radiation 

An argument often made against the beaming model of PULXs is the belief that the observed high (< 45%) pulsed 
fraction of the X-ray emission cannot appear in radiation passing through a funnel because reflections from the funnel 
walls would destroy the signal coherence. 

We note first that although many (if not most) ULXs probably contain magnetised neutron stars, only a small frac- 
tion show pulses. Second, (2020) point out that the argument above assumes (without any explicit 
statement) that the spin axis of the neutron star is aligned with the accretion disc axis, and so with the beam anid In 
fact observations of high—mass X-ray binaries show that the neutron star spin is generally not aligned with the disc 
axis, probably because of the asymmetric nature of the supernova giving birth to the neutron star. 

The beamed emission will appear sinusoidal, (i.e. the light curve has no constant intervals) if three conditions hold 
simultaneously: 

(a) there is significant misalignment of spin and beam (disc) axes (cf Figure[i5), and 

(b) the emitting region is not essentially symmetrical around the spin axis, and 

(c) its linear scale is comparable to the neutron star radius, so that a significant fraction of its area, but not the whole 
of it, is occluded by this star as it spins 


i ~ 45° - 90° 
m ~ 0° - 45° 


i~ 0° - 45° 
m ~ 45° - 90° 


pulses seen m BO (weak) pulses 


Figure 15 Effect of spin orientation on pulsing. Left: the neutron—star spin and the accretion disc beaming axes are 
strongly misaligned, so that a significant part of the pulsed emission is periodically occulted by the neutron—star body 
at certain spin phases. This gives a large pulse fraction. Right: the neutron-star spin and central disc axes are assumed 
to be substantially aligned (this special condition is assumed without any explicit statement in papers asserting that 
beaming cannot produce pulsed emission). Much of the primary X-ray emission is scattered by the walls of the 
beaming ‘funnel’ before escaping. In this kind of configuration the pulse fraction is reduced or suppressed. 


2020). 


10This unstated assumption has a similar status to the implicit belief that ULX nebulae must be Strömgren spheres powered by the direct radiation 
of the ULX, rather than wind bubbles made by its powerful near—spherical wind — see Sections and In both cases the effect is to ‘rule 
out’ beaming. 
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Two processes can align the spin and disc (beam) axes. One is warping of the central disc into the neutron-star 
spin plane through differential torques (both magnetic and precessional), and the other is direct accretion of angular 
momentum (characterized by Ry) from an unwarped disc lying in the binary orbital plane. The latter seems more 
likely, as any interruption in the accretion flow means that disc warping has to ‘start again’. Neither process will be 
be very efficient in the Be-star PULXs, where accretion rate is less than < 10Mgaq and confined to short transient 
episodes. 

The outcomes of the two spin orientations are very different. 


e If there is strong misalignment of the spin and beaming axes, a significant part of the pulsed emission can 
escape without scattering, giving a large pulse fraction (Fig. left). This must be maximal at the highest 
X-ray energies, as scattering makes the X-rays both softer and less pulsed. This correlation of pulse fraction 


and X-ray energy is well known for observed PULXs 2017). 


If instead there is substantial alignment of the spin and central disc axes, much of the primary X-ray emission 
is scattered by the walls of the beaming ‘funnel’ before escaping (Fig. right). Since the light-travel time 
across the funnel is usually comparable with the pulse duration, the pulse amplitude can be severely reduced. 
This would result in an inability to detect many highly collimated PULXs through pulse searches 
fet _al.|/2019); as accumulate and improve in quality, this picture will be tested more thoroughly. The pulse 
fraction can of course also be reduced or entirely removed if enough matter accretes to weaken the neutron—star 
magnetic field. This presumably happened in the case of Cyg X—2, which is a survivor of a phase of strongly 


super—Eddington accretion (King and Ritter||1999). Here the neutron star is not noticeably magnetic, and has 
probably gained a few tenths of Mo during the super—Eddington phase. 


These outcomes appear to be consistent with observations of ULXs. Very few ULXs show pulsing. But if most 
ULXs are the direct descendants of high—mass X-ray binaries (HMXBs, once the companion star fills its Roche lobe), 
or Be—star HMXBs, the vast (unpulsed) majority must contain neutron stars. First, almost all HMXBs contain neutron 
stars rather than black holes (2017). Second, neutron-star systems are more 
super—Eddington and more beamed than black—hole systems with the same mass transfer rate, and therefore have 
higher apparent ULX luminosities (see Section 3.6). In addition, there is at least one ULX whose spectrum shows a 
CRSF corresponding to a pulsar-strength magnetic field (Middleton et al.|/2021b), but which does not show pulsing, 
strongly suggesting that is it an aligned accretor. 

The vast majority of ULXs do not have detected pulses, despite containing neutron stars. Unless this is caused by 
a lack of data for ULX pulse searches it implies that alignment of spin and central disc axes is rapid, and possibly that 
field suppression through accretion occurs also. As expected in this picture, Be—star systems are prominent among 
the PULXs since accretion rates are relatively low and the episodes of increased luminosity are transient. 


3.8. Numerical simulation of super-Eddington accretion flows 


Most of the numerical simulations of super-Eddington accretion on to compact bodies concern black holes. The 
purely absorbing inner boundary conditions make the problem simpler than for neutron stars, which have hard surfaces 
and, in the most interesting case, strong and rotating magnetic fields. We first review the simulations of black—hole 
super-Eddington accretion. 


3.8.1. Numerical simulations of black-hole accretion 

Analytical or quasi-analytical solutions of super—Eddington accretion on to black holes represent two extreme 
possibilities: either the local emission is kept close to the Eddington value by blowing out the “excess” matter (Shakura| 
1973), or by advecting the “excess” energy towards the black hole, where it disappears from the Universe 
(Abramowicz et al.|{1988). In both cases the total luminosity radiated by the accretion flow is 1 + Insio times the 
Eddington luminosity, where mo corresponds either to a constant accretion rate (slim disc) or is equal to the rate at 
which matter is brought to the accretion flow, the accretion rate decreasing to keep the local flux at the Eddington 
value (wind/outflow solution). Although the accretion rate reaching the black hole can differ between the two cases 
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by several orders of magnitude, in neither case can the emitted luminosity be larger than ~ 10 Lead!" Then for a ULX 
with Ly > 10 Lgaa, this apparent luminosity cannot be equal to the physical luminosity, whatever the accretion flow 
model: the observed super-Eddington luminosity must be beamed. Both extreme models imply beaming depending 
on the accretion rate. 


For slim discs, (Wiktorowicz et al.||2017) fitted the relative semi-height of the scattering photosphere calculated 
analytically and numerically by (2016) by the formula: 


pho 1.6 
ee, (87) 
R 1+4 
and from i 
Aphot \ 0 
(=) = tan 5, (88) 
obtained the beaming factor 
0 
b = 1 — cos 5" (89) 


This beaming becomes nearly constant for m > 100 which gives the maximum beaming for b ~ 0.15. For the 


~N 


windy |Shakura and Sunyaev| (1973) model with beaming given by Eq. (55). b = 0.15 corresponds to m = 22. 


One can expect that in reality the solution describing super-Eddington accretion on to a black hole may have both 
outflows and advection, but the relative proportions can only emerge from numerical calculations and comparison 
with observations. Unfortunately both these methods encounter serious difficulties. 

Before we discuss numerical solutions describing super-Eddington mass supply to black holes, we should note 
that, without estimates of the black-hole masses in ULXs, these models remain untested, since we do not know the 
true Eddington luminosities. We now know that the masses of black holes descending from stellar evolution can differ 
by an order of magnitude (Abbott et al.|/2020), so this is a real handicap. Currently there is only one upper limit on 
the compact-object mass in a ULX (Motch et al.| [2014), and here the accretor was later found to be a neutron star 
[Fürst et al.|/2016). 

In discussing numerical simulations of super-Eddington accretion flows, we emphasize that these cannot so far 
represent the conditions expected in ULXs, i.e. a disc of at least > 10° R, in size, geometrically thin for R = 15mRg. 
This is not a surprising result: MHD codes attempting to model the disc viscosity from a first-principles treatment of 
the magnetorotational instability currently find values of the dimensionless viscosity parameter a significantly smaller 
(< 0.01) than observational estimates from time variations, which strongly suggest a = 0.3 (Tetarenko et al.|/2018). 
The reason for this discrepancy is likely to be numerical diffusivity, as already recognized in the solar and MHD-fluids 
literature. One should add, however, that the viscosity parameter can be reliably determined in the outer disc’s regions 
only (see, e-..Kotko and Lasotal2012}. 

Most, if not all, numerical simulations of such flows start with a torus, with an assumed distribution of angular 
momentum that after relaxation forms an accretion flow. In most cases the “circularization radius” of a torus—generated 
flow (the Keplerian orbit corresponding to the torus angular momentum) is kept rather small (~ 30 — 50 R,). This is 
because the simulations aim to get quasi-stationary solutions, but this requires run times significantly longer than the 
viscous time of the accretion flow. 

Since the viscous time at radius R is 


_/H —2 RÈ 1/2 
fis ~ a"! (5) (z , (90) 
one can write 
s| R Vera \!(H/RY? 
vis ~ 1.1x1 — — Ío; 1 
f eee (car) (a7) (5%) s CA 


where the dynamical time t, = R,/c, is often used to express the duration of numerical calculations (t, = 5 x 1075 
s, for a 10 Mo black hole). Typical 3D GRRMHD (general-relativistic radiation magneto-hydrodynamic) calculations 


‘For “hyper-Eddington” accretion on to an intermediate-mass or super-massive black hole with ri = 5000, the luminosity is still ~ 30 Lgaa 


(Sakurai et al. [2016). 
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by |Sadowski and Narayan} (2016) require runs of 20000 fz, so it is not surprising that their quasi-stationary solutions 


require rather small discs, with R < 40Rg. 

The pioneering simulations of super-Eddington accretion by[Eggum et al.](1988) were able to represent only 0.6 
seconds of physical evolution of the system, while the improved scheme of (2002) increased this duration 
to 1.6 seconds, both shorter than the corresponding viscous time. 2005), using 2D radiation hydro- 
dynamic (RHD) simulations were the first to obtain a quasi-steady structure of a supercritical accretion flow onto a 
black hole. It took another 10 years for other groups to join the Japanese team of Ken Ohsuga in trying to model 


supercritical accretion onto compact bodies. |Sadowski and Narayan] (2015); Sadowski and Narayan) (2016) used a 
general-relativistic RMHD code first in two, then in three dimensions, while|Jiang et al.|(2014) performed 3D simu- 


lations in the framework of the so—called pseudo—Newtonian potential of/Paczynsky and Wiita 1980)]! Except for 
Takahashi et al.|/(2016), who used the general-relativistic description of the black-hole accretor, other simulations of 


the Ohsuga group (Ohsuga et al.|/2005} |Ohsuga and Mineshige} |2011 2015 2018 


2021) used the pseudo-Newtonian description of particle motion in the gravitational field of the accreting body. All 


Table 5 Numerical Models of super-Eddington black hole accretion (Adapted from 2021) 


Paper Method Compton | Rout Rx Rass Rap Megu Moutfiow 
[Yes/No] | [R] | [Rel | [R] | (Rel |b] |i 
Kitaki+21 2D-RHD Yes 6000 4860 | ~ 1200 | ~ 540 ~ 18 ~ 2.4 
Ohsuga+05 2D-RHD No 1000 200 ~ 60 ~ 400 ~ 260 
Ohsuga+11 2D-RMHD No 105 40 ~ 10 ~ 150 ~ 100 
Jiang+14 3D-RMHD No 100 50 ~ 40 ~ 660 ~ 22 ~ 40 
Sadowski+15 2D-GR-RMHD | Yes 5000 42 ~ 70 ~ 1280 | ~ 42 ~ 700 
Sadowski+16 3D-GR-RMHD | Yes 1000 40 ~ 20 ~ 520 ~ 18 ~ 52 
Hashizume+15 | 2D-RHD No 10000 | 200 ~ 200 ~ 460 ~ 15 ~ 50 
Takahashi+ 16 3D-GR-RMHD | No 250 34 ~ 20 ~ 600 ~ 20 
Kitaki+18 2D-RHD Yes 6000 600 ~ 400 ~ 840 ~ 28 ~ 30 
Jiang+19* 3D-RMHD Yes 1600 80 ~ 30 ~ 760 ~ 25 


Rout — radius at the outer boundary, Rx — initial Keplerian radius (“circularization radius”), Rgss — radius, inside which a quasi-steady state 
is established, Rtrap — photon-trapping radius derived based on equation|34)(with right hand side multiplied by ~1.5), Meu is the accretion 
rate onto the black hole, Moutfow is the outflow rate at around Rout. It is indicated whether the Compton scattering effect is taken into 
account or not. 


* In these simulations the black-hole mass is 5 x 108 Mo 


these simulations produce rather similar results: strong outflows and regions of the inflow dominated by advection. 
They differ mostly in details, except for|Jiang et al.|(2014) who observe an effect not seen by the other teams: vertical 
advection of radiation caused by magnetic buoyancy that transports energy faster than photon diffusion. This allows 
a significant fraction of the photons to escape from the surface of the disk before being advected into the black hole; 
vertical “advection” reduces horizontal advection. The reason for this discrepancy between simulation results is not 
understood. The main suspect is the treatment of radiation in the code. While|Jiang et al.|(2014) solve the full radiative 
diffusion equations, other authors use approximations, such as flux limited diffusion (FLD) method or the M1 closure. 
This might explain differences between the radiative energy distribution and radiation collimation observed between 
the two sets of simulations. On the other hand, in (2014) the inner edge of the simulation box is located 
outside the horizon which, as showed A eS ee oan can lead to reflection of energy that would 
have normally been advected into the black hole. 


However, as pointed out and discussed in some detail by |Kitaki et al.| (2021), all the simulations above suffer 
from one serious drawback that puts into question their relevance to direct modelling of ULXs: their circularization 


!2Since the Paczyriski-Wiita “potential” mimics the general-relativistic description of massive particle orbits in the black—hole’s equatorial 
plane, it should actually be called “pseudo-telativistic’’, but it is too late to change the widely accepted nomenclature. (The name “Paczyński” was 
misspelled in the original publication, hence the spelling used in the reference list.) 
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Figure 16 Schematic view of the structure of the super-Eddington accretion flow and associated outflow based on our 


numerical results. The black arrows indicate the gas motion. (Kitaki et al.||2021) 


radii, as well as the outer radii of the quasi-steady disc region are always inside the trapping and spherization radii 
(see Table[5p. Therefore although all numerical models predict wind plus advection dominated flow configurations it 
remains to be seen where such features occur, and what their properties are in real systems with, say, Rout > 10*R, 
and Rsph > 750R, (corresponding to Rout > 3 X 10!° cm for a 10 Mo black hole and mn > 10, respectively). 

have made a substantial contribution towards reaching this objective (see first line of Table[5). 
They performed simulations with an outer radius of 6000R, and circularization radius 4860R,. Their calculations 
achieved a quasi-stationary state inside Rąqs ~ 1200R,, a radius 3 to 120 times larger than in previous simulations of 
super-Eddington accretion flows. used a 2D R-RH code, with a viscosity parameter œ = 0.1. The 
black hole mass is 10 Mg and a Paczyriski-Wiita “potential” is used to describe massive—particle motions. Matter is 
continuously added at Rou at a rate ń = 70, with specific angular momentum 3.12 x 10'8cm?s~!. All matter arriving 
at Rin = 4R; is absorbed. 

The structure obtained by|Kitaki et al.](2021) is schematically shown in Figure [16] The accretion rate onto the 
black hole is mpgy = 18. The flow forms a well—delimited disc structure with a surface defined as the loci at which the 
radiation force balances the gravity in the radial direction. Near the black hole and up to R ~ Rrap the disc height H is 
proportional to H (H/R ~ 1), above this radius it is constant H ~ 22—347ngyR,. The trapping radius Rrap ~ SOmpyRg. 
This means, as expected from simple models (Section[3.2), that for R S Rap the flow is advection dominated, while at 
larger radii it behaves like a radiation—pressure dominated, Shakura-Sunyaev-—type accretion disc. Such discs should 
be thermally unstable, but this does not seem to occur in|Kitaki et al.](2021) simulations despite the calculation time 
amply allowing the instability to develop. The reasons for this remain unclear (see, however, 2013) but 
are worth investigating. 

Closer to the black hole, at Riau ~ 26mguRg, an outflow is produced, but the expelled gas is not able to reach 
infinity and falls back on to the disc. Only at RẸ, ~ I5ripyR, is a real outflow (wind) blown out. The formula for 
Rọ, has the same form as that for Rspn (Eq. but with py instead of 7 which in Eq. corresponds to the rate 
at which matter is fed to the outer radius (in the simulations under discussion, m = 70, while mgy = 18). We discuss 
this difference later in this section. Finally, in the regions closest to the black hole, below Rinfow ~ 4.47ipHR,, most of 
the gas flows in at a rate sn = 18, with some weak outflow now, ~ 1.3. 
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The|Kitaki et al.|(2021) solution is fully consistent with what one expects from an accretion-launched wind. The 
wind speed is v ~ 0.1c, and Lmech ~ 0.05 — 0.08Lraq, with Drag = Ly ~ 2 — 3Legaa, SO Lmech © 0.1 — 0.24Lgaa. The 
launch radius of the wind is Riau ~ few 100R,, in agreement with the wind speed v ~ (2GM/Ria)'/?. This gives a 
wind momentum rate Mout = (2/v) Lmech ~ 20Lmech/C ~ Lrag/¢ as one would expect when the mass supply is not 
hyper-Eddington (here it is ~ 7 = 18). These values are similar to analytic estimates 2003) for 
momentum—driven winds. 

This solution suggests that for black-hole accretors, below R ~ 1000R,, most of the mildly (m ~ few x 10) 
super-Eddington mass accretion is simply swallowed, as predicted by slim-disc models. The outflow to accretion ratio 
Mitouttlow/MMgu — 0.14 is much lower than in previous simulations which |Kitaki et al.|(2021) attribute to the sizes of the 
simulated flows, theirs being much larger than in the preceding calculations, thus avoiding having the whole simulated 
domain in the region where the flow is puffed—up from the start. But we note that although Rqss ~ 1200R, is 20 to 
30 times larger than the corresponding radius in R-GR-MHD calculations, it is still “only” 1.8 x 10°cm, well inside 
the circularization radius. The net accretion flow for R < Rgss is roughly constant (m ~ 18) while the mass inflow 
rate decreases from m ~ 30 at R ~ Rg; tom = 18 for R ~ 40R, (m ~ R-°), but it is not clear what happens above 
1200R,, where the simulated configuration is not relaxed. 

For example, for R > 2012R,, the gas pressure dominates over radiation pressure and the disc should be geo- 
metrically thin (see Eq. with m = 70). This should reduce the outflows, thus making the whole configuration 
potentially inconsistent. It seems that we are still some distance from having a complete description of the whole 
picture of super-Eddington accretion onto black holes. 

Indeed, recently, simulated super-Eddington accretion onto a 10° Mg black hole. The physical 
outer radius of their simulation domain is R = 1500R,, four times shorter then in (2021). In their case, 
however, the steady-state radius Rgs; = 1500Rg, slightly longer than in{Kitaki et al.| (2018). The flows they study settle 
down to a quasi-steady state in millions of the orbital timescale, which in their case corresponds to the viscous time at 


~ 1000R, (they take œ = 0.01). In their case Min œ r” with p ~ 0.5 — 0.7 and contrary to the (2021) only 


a small fraction of the inflowing matter feeds the central black hole. Indeed, their solutions can be represented as 


Meu = 17.1M (= F 
BH = . Edd 300 
Poa = 0.51cMeaa (2%) (92) 
out = U.ILC vaa (Z) 
: 0.5 
mo 
2 175(22) -1 
B 300 
and the wind velocity is 
17x 103e( Lt.) 93 
wind) = 1./ X aA > 
(wind ) (2) (93) 


where 8B = Mow /Min is the wind loading fraction. It is only ~ 0.13 in (2021) but ~ 20 in the model of|Hu 
(2022), so that even in the case of a black hole, advection does not play a significant role when accretion rate 


are very strongly super-Eddington. Therefore the main difference between the two numerical models is that in|Kitaki] 
(2021) most of the outflowing gas fails to escape from the outer boundary and falls back onto the disc, while 
in most matter is lost in the wind. (2022) attribute the large differences between the two 
approaches to several factors, numerical and physical. First, impose in their computational domain 
equatorial mirror-symmetry of the accretion flow across the equatorial plane. It seems that this imposed symmetry 
tends to suppress global convective motion that crosses the equator and is known to produce outflows coherently 
through the equator rather than the polar eT EDO] 2013). The second reason would be that {Kitaki et al.| 
assume a viscous parameter of œ = 0.1 (Hu et al. assume «œ = 0.01) which is apparently known to cease 
convective/turbulent motion of advection dominated flows and thus to reduce the outflow rate (see, e.g. {Inayoshi et al.} 
2018). These differences do not necessarily favour the model with stronger outflows. For example, it seems that high 
values of the viscosity parameter œ > 0.1 are appropriate to the description of high—accretion rate flows 
fet al.| [2018). Only a self-consistent calculation with a MRI-determined (anomal) viscosity could decide if numerical 
models can reproduce the mechanical power observed in some ULXs. In the most recent paper of the Ohsuga 


group (Yoshioka et al.}|2022) it is concluded that the fraction of the failed outflow decreases with the decrease of 
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Min. They leave the comparison with the|Hu et al.|(2022) simulation for a future paper. What’s more important, their 
simulations confirm the effect of radiation beaming at high accretion rates since they found that “the higher Min is, the 
more vertically inflated becomes the disk surface, which makes radiation fields more confined in the region around 
the rotation axis”. 

An interesting result of most numerical simulations is that they do not see any evidence of the photon bubbles 
proposed as a solution of the super-Eddington luminosity problem by (2002). When such bubbles do 
appear, they play no significant role in super-Eddington accretion flows (see, e.g., (2014). 

We are even further away from understanding super-critical accretion through numerical simulation when the 
accretor is a rotating, magnetized neutron star, as discussed in the next section. 


3.8.2. Numerical simulations relevant to PULXs 

The difficulty of simulating super—Eddington accretion on to rotating, magnetized neutron stars is probably best 
illustrated by the fact that until now only four papers have dealt with the problem: [Takahashi and Ohsugal (2017); 
[Takahashi etal (2018}[Abarca etal 2018] 2021). The maximum dipole field strength considered is ~ 10°G (in two 
cases B = 0 G), and none of the simulations describes a rotating accretor with a hard surface. 

All attempts to describe super—Eddington accretion on to neutron stars use GRRMHD codes, similar to, or inspired 
by, the codes used in the description of black—hole accretion, described above. There are, however, two further 


difficulties that must be overcome when such codes are applied to magnetized neutron stars. First, it is not easy to 
describe in the same code both magnetically—dominated flows with a force-free character and the “normal” properties 


expected from the (GRR)MHD equations. This difficulty has been partially removed by |Parfrey and Tchekhovskoy 
(2017). 


The second difficulty is in fixing the inner boundary conditions, i.e. the physical properties of the accreting matter 
at the the neutron star surface. In contrast to black—holes, whose surface is always fully absorbing independent of 
what — and how much — is falling onto it, the properties of matter arriving at the neutron-star surface depend on the 
magnetic field strength and the accretion rate (see|Abarca et al.|/2021| and references therein). As in most blackhole 
simulations, limits on the available computational time restrict the range of radii for which the solutions relax to quasi— 
stationarity; in all simulations Rass < 50R,, while the spherization radius is Rsph << 500R, and the computational 
domain extends up to 1000R,. 


Takahashi et al.|(2018) and (2018) consider the case of super—-Eddington (m ~ 10) accretion onto a 


non-magnetic and non-rotating neutron star with a reflective surface, and compare it with the case of accretion onto a 
black hole. They both observe more powerful outflows in the neutron-star case and {Takahashi et al.](2018) find that 
for a neutron-star accretor, the accretion flow inside R < 10R, (in their case Ross ~ 10R,) corresponds to the SS73 
windy solution (Sect. [3.3}), whereas in the black-hole case, the flow is advection dominated!*} The strong wind blown 
out from the vicinity of a super—Eddington accreting neutron star is very optically thick to electron scattering “which 
would lead to the obscuring of any NS pulsations observed in corresponding ultraluminous X-ray sources” 
fet al.| [2018). Since the model of a non-rotating and unmagnetized neutron star would produce no pulses in any case, 
this is not a big drawback but stresses that the presence of optically—thick outflows from super—Eddington accreting 
neutron stars cannot be neglected when describing their X—ray emission. 

performed a 2D axisymmetric radiative GRMHD simulation of accretion onto a neutron star 
with a 2 x 10'° G dipolar magnetic field. The combination of the hard surface and confinement of the gas into an 
accretion column near the star, allows the flow to release radiative energy at a rate of several times the Eddington 
limit. They compared their results to the KLK17 model and found a larger beaming intensity, but they believe that 
post-processing would show a less intensely beamed distribution of radiation at infinity. They also note that they do 
not model the same system as considered in KLK17 (e.g. they have super-Eddinton accretion rates in the accretion 
column). In their case the distance between the Alfvén radius and spherization radius is large, and their star is non- 
rotating. This paper reveals the difficulties current numerical simulations still face in modelling real PULXs. 

In all simulations, the estimated minimum values of the inverse of the beaming parameter b (maximum beaming) 
obtained from various models range from 1/b = 1/30 — 1/20 =% 0.033 — 0.05, which compares quite well with the 
range of values b = 0.009 — 0.6 from Table[3]in Section [3.7.1] i.e. with the results of the KLK17 model. 


This clear physical difference between accretion flows onto compact objects which have, or do not have, surfaces, makes the advection- 


dominated solutions of|Chashkina et al.|(2019) rather unrealistic. 
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Figure 17 Stream lines around a (1.4 Mo) neutron star (left) and a (10 Mo) black hole (right). Red (blue) lines indicate 
that the radial velocity is in the positive (negative) direction. (From|Takahashi et al.|/2018). 
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The simulations of super—Eddington accretion onto neutron stars thus provide strong support for the assumptions 
of the KLK17 PULX model, i.e. that the accretion flow inside the spherization radius is well described by the SS73 
windy model, and that the X-ray luminosity is beamed. Although these simulations assume weakly magnetized and 
non-rotating neutron stars, their basic results: weak advection, unlike in the black—hole case, and therefore strong 
optically thick outflows, probably hold for stronger fields and faster rotation of the accretor. 


3.8.3. Numerical simulations of PULX accretion columns 

Numerical models of super-Eddington accretion onto a strongly magnetized neutron star are given by [Kawashimal| 
hereafter K16) and{Kawashima and Ohsuga|(2020] hereafter K20). The simulations describe the last part 
of the accretion flow when the infalling plasma forms a column above the neutron star surface. In both papers, the 
accretion column is represented by a truncated cone with a base (“polecap”) occupying a fraction 0.07 of the stellar 
surface. Since the method used is radiation hydrodynamics, the value of the magnetic field does not appear explicitly, 
but in K16 it is estimated to be ~ 10!°G, while in K20 the magnetic field strength is > 10!?G. In the first case, the 
plasma is allowed to move transversely, in the second, the motion of matter is restricted to the radial directions. In 
K16 the initial uniform density in the column gives a time-averaged accretion rate of m = 50; K20 uses the same 
initial set-up but also a lower-density model, giving 7 = 3. In both types of model, a radiative shock forms above 
the neutron-star surface (fixed at 10 km) at a height ~ 3 km for the ”weak magnetic field” and at ~ 10km in the case 
where the accreting matter motion is restricted to the radial direction. In the “weak-field” simulations, practically 
all the radiation is released below the shock and radiated sideways. The apparent luminosity is highly anisotropic; 
it can reach ~ 50 Lega sideways, but is close to Eddington vertically along the accretion column. In contrast to the 
“weak-field” case, where the radiation field below the shock is homogenized by circular motions of radiation bubbles, 
the “strong-field” simulations show that matter is accreted along the side walls of “hollow cones” while matter inside 
the cones is blown away by radiation. For 7 = 3, only one hollow cone forms, while for very high accretion rates 
(m = 50), the column contains three such structures. Also in this case the apparent luminosity is highly anisotropic, 
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Figure 18 Schematic picture explaining energy flow from gas (potential energy) to outgoing radiation within a super- 
critical accretion column. The blue arrows represent the energy flow carried by gas: their length and width are drawn 
in proportion to the kinetic energy flux and the mass accretion rates, respectively, whereas the red arrows represent 
energy flow carried by radiation and their widths are drawn in proportion to the radiation energy flux 


radiation still managing to get out mainly sideways, but radiated from a larger surface than in the “weak-field” case. 
The sideways luminosity can reach 30 Lgaa. 

Using the K16 model, {Inoue et al.|(2020) calculated the light-curves produced by such an accretion column, taking 
into account general-relativistic effects on light propagation near the neutron-star surface. The pulses they obtain are 
quasi-sinusoidal. They obtain a high pulsed fraction of the emitted light, from 5% up to 50%, depending on the 
magnetic and observational inclination angles. However, a narrow magnetic column for a highly super-Eddington 
accretion rate, and field of ~ 10!°G does not seem to give the best representation of a real PULX since high accretion 
rates and low magnetic fields should instead favour wide accretion columns. 

The papers K16 and K20 are important contributions to the understanding of supercritical accretion in strong 
magnetic fields, but it is unlikely that they correspond to the accretion structure of real PULXs. The accretion flow 
in these systems probably does not form a column extending to 2 x 108 cm (~ 200 neutron-star radii) above the 
neutron star surface. Since PULXs are binary systems, matter transferred from the neutron-star’s companion forms 
a disc extending down to a magnetosphere where the pressure of the accreting matter is comparable to the magnetic 
pressure. The size of the magnetosphere will determine the angular size of the column, which will vary with magnetic— 
field intensity, differing by two orders of magnitude. For the parameters of K16, the magnetospheric radius is close to 
the neutron star surface at ~ 3 x 10° cm (see Eq. and for the K20 set-up Ry < 108cm. What is most important, 
however, is that accretion within the magnetosphere is most likely highly dissipative, since the flow along field-lines 
requires a highly supersonic flow to bend and shock, while in K16 and K20, the first shock occurs very near the 
neutron-star surface, because in the simulations (different from the schematic view of Fig. the field-lines are 
straight. By necessity, the K16 and K20 simulate only the last part of the PULX accretion flow and they assume that 


it arrives at the magnetosphere at a super-Eddington rate, which is rather unlikely in reality (Takahashi et al.}|2018). 


3.9. PULXs as magnetars in binary systems 


The discovery that ULX M82 X-2, with an apparent luminosity L = 1.8 x 10*° erg s™!, is a pulsar with a coherent 
periodicity P = 1.37 s (Bachetti et al.||2014), raised the possibility that its apparent super—Eddington luminosity might 
be intrinsic (cf [Tong] 2015} [2015a). If the strength of 
its magnetic field is magnetar-like (B ~ 10'*G), the reduction in electron scattering opacity for some directions and 
polarizations can be quite substantial. 


At first sight, this idea is similar to the one motivating (1992) to explain the luminosity L ~ 10* Leasa 
observed in the soft-gamma repeater (SGR), GB 790305 (“the 5 March event”), by the presence of a magnetic field B = 
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3 x 10'4G. Ina very strong magnetic field, the Thomson and Compton scattering cross-sections are strongly reduced 
for photons with energies E, much lower than the energy corresponding to cyclotron frequency Eeye = 11.6B12 keV, 


where B1210! G = B (1971). From[Herola|{1979}) one has 


T1 2 E, ? 2 
— xsin“ + (=) cos’ 0 (94) 
OT cyc 

E,\ E 
Za z) , for Y x 1, (95) 
OT cyc cyc 


where indices 1 and 2 correspond to the two linear photon polarisations, and 6 is the angle between the directions of 


the magnetic field and light propagation. (1992) gives an approximate relation between the “new” critical 
luminosity and the Eddington luminosity, defined through the Rosseland mean opacity: 


Lorit 2 | E jo 


Lea 2 \2 x 10!4em s72 Oe) 


where g = GM/R, when Leit >> Lega. This critical luminosity is of course only approximate because it is not 
properly integrated over angles, but estimates of the Rosseland mean opacities vary with angle only by a factor two 


(Paczynski}| 1992). 
From Eq. (96), for a 1.4 Mo neutron star, a sub-critical luminosity L < 0.1 Lerit implies 


4/3 
Va 2% 10"( ergs /, (97) 


10¢G 
so that explaining PULXs with L > 10*° ergs! assuming isotropic radiation, requires magnetic field strengths that 
would be considered extreme, even among magnetars. 

In PULX models by [Tong] (2015); Dall" Osso et a (2015): [Ekşi etal 2013) OTT), the fields 
are in the range ~ 3 x 10? — 6x 10°° G (the last value considered by the authors to be “rather unphysical”), so that the 
critical luminosities corresponding to the first three PULXs discovered put them safely into the sub-critical luminosity 
regime, even with a magnetic moment derived by magnetic torque considerations. 

But there is a fundamental difference between SGRs and PULXs: SGRs are powered by magnetic energy 
[and Thompson] [1992}, while the luminosity of PULXs comes from accretion. For SGRs the magnetic field allowing 
L = 10° Lgaa is consistent with pulsar slowdown and the required magnetic energy. There is no known alternative to 
magnetar-strength fields in SGRs and anomalous X-ray pulsars, and there is so far no observational evidence against 
the presence of fields with these strengths. 

This is not true of the magnetar hypothesis for PULXs. First, no magnetar is observed to be a member of a binary 


system. All ~ 30 known magnetars (or “candidate” magnetars) are isolated neutron stars (Olausen and Kaspi}|2014). 


Conversely, all neutron stars in X-ray binaries are observed to have magnetic field strengths < 10'°G, below the 


magnetar range (see, e.g.,/Revnivtsev and Mereghetti} |2016). 


The super-strong magnetic field in magnetars is thought to be formed in a process which destroys (or merges) the 
binary, of which its progenitor was a member 2016). This would explain their absence in binary systems. 


Magnetars are supposed to form through neutron-star mergers (e.g., 2013} 2017) 


or in type Ibc supernova outbursts in massive binaries. The first mechanism obviously produces a single object. In the 
second, the superluminous supernovae disrupt the binary and produce an isolated magnetar, as demonstrated by the 


magnetar CXOU J1647-45 in the young massive cluster, Westerlund 1 (Clark et al.|/2014). 
We note that “most of theoretical considerations do not favour even existence, not speaking about active decay, of 


magnetar-scale fields in neutron stars older than 10° yrs” 2022). 


These considerations make the presence of rapidly rotating, strongly accreting magnetars in binary systems highly 
unlikely. Their presence in ULXs would evidently require a cosmic conspiracy. 


3.10. ULX populations 


The beaming formula allows one to study the populations of disc-wind beamed ULXs in galaxies. If the host 
galaxies have space density n, Mpc~3, and each host contains N ULXs with randomly oriented beams, one would 
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need to search through ~ 1/Nb galaxies, occupying a space volume ~ 1/n,Nb, to be sure of being in the beam of one 
ULX. Then the nearest ULX must be at a distance 


D ~ OTN)! n? Mpe (98) 


1/3 
ans (s) 
ń where m = m/10, and has apparent (isotropic) luminosity 
L=22X 10m, m? erg s7! (99) 


where m, = m/(10 Mo). The scalings m,, mı are appropriate for black holes or mass ~ 10 Mo, so by the argument 
following Eq. the luminosity normalization here would be 2.2 x 10*° ergs™! for neutron star ULXs. These 
relations can be used to show that the ULX luminosity function of the Local Group is consistent with the idea that 
beamed ULXs are the evolutionary stage following the standard high-mass X-ray binary phase (cf 


OTO) 


3.10.1. Pseudoblazars? 
The beaming described by potentially allows one to see suitably—oriented ULXs out to very large distances. 


The extreme object SS433 has sn ~ 300 — 104 (King et al.}|2000b}|Begelman et al.|/2006| see Section|3.14.1]). With 


tn = 10*ring the nearest object of this kind which we would observe within the beam is at a distance 
3 1/3 
Dmin ~ (= | ~ 660N7"7 rn? M 100 
Fane ;] Ma ps (100) 


where we have taken ng ~ 0.02 Mpc™ as appropriate for L* galaxies. The apparent isotropic luminosity of such an 
object would be 
Lspn = 2.2 x 10° mm, ergs™". (101) 


This distance and apparent luminosity are typical for active galactic nuclei (AGN), so without further information one 
could not easily pick out a bright ULX like this from a sample of observed AGN. But unlike a genuine AGN, there is 
no reason to suppose that a given ULX lies in the nucleus of the host galaxy. There are possible candidates for such 
systems: the BL Lac system PKS 1413+135 2002) has distance D ~ 1 Gpe and isotropic luminosity 
~ 10“ ergs™! , but lies at 13 + 4 mas from the centre of the host galaxy, and HLX-1 could be a system of this type 
(cf Section 2.7). Of course it is difficult to exclude the possibility that this system may lie at the nucleus of a fainter 
line-of-sight galaxy. 


3.11, ULXs and binary stellar evolution 


Sections [3.6] [3.7]and[3.8]above have shown how the phenomena characterizing ULXs arise, in many cases, from 
strongly super—Eddington mass transfer in binaries containing compact objects (white dwarfs, neutron stars or black 
holes) which have reached the end of their stellar evolution. We should now ask what causes these very high mass 
transfer rates. To do this we need to consider the evolution of binary systems where one of the stars is compact (a 
black hole, neutron star, or white dwarf) and the companion star transfers mass to it, powering accretion. The rate of 
mass transfer is determined by the evolution, both of the companion star, and of the binary orbit, as mass transfer in 
general changes this. Dissipation in the mass transfer process generally makes the binary orbit circularize. 

Kepler’s laws give the orbital angular momentum of a circular binary system consisting of stars of mass Mı, M2 
as 


1/2 
Sa) l (102) 


J=M\M) (= 
where M = Mı + M3 is the total binary mass, and a is the orbital separation. We take star 1 as the compact object, and 
star 2 as an extended star. In a high—-mass X-ray binary (HMXB), star 2 is hot, and loses mass in a wind. The X-rays 
result when some of this wind is captured by the compact star (a neutron star or black hole). But as the extended star 
2 approaches the end of core hydrogen—burning it begins to expand. Eventually its radius R fills its Roche lobe 


Rz = f(q)a (103) 
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(where q = M>/M, is the binary mass ratio, and f(g) < 1 is a slowly—varying function of q). The star begins to 
transfer mass through the inner Lagrange point, the saddle point of the combined gravitational-centrifugal, Roche 
pseudopotential between the two stars. It is generally a good approximation to assume that all the mass lost by star 
2 is accreted by the compact star 1, and that the orbital angular momentum J is conserved in this process. Then 
J = M =0,M, + M, = M = 0, and logarithmic differentiation of Eq. gives 


i= Al (104) 


Since -M > 0, we see that if the extended, mass-—losing star is more massive than the compact star, the effect of 
mass transfer is to shrink the binary separation (å < 0). As f(q) varies only weakly with mass ratio q, Rz also shrinks, 
tending to increase the rate of mass transfer. The result is that mass transfer from the more massive to the less massive 
star is self-sustaining, at least for a time. 

But there is a limit to this process. The outermost stellar gas — above the saddle point of the Roche potential — 
is lost on a dynamical timescale, and gas from below flows towards the saddle point on the same timescale, trying 
to refill the Roche lobe. But this gas is not now in thermal equilbrium because its rapid rise means that it retains the 
same specific entropy it had when deeper in the star. But in stars with a radiative envelope (typically, blue stars with 
effective temperatures Ter > 104 K), the equilibrium entropy increases outwards, so the new gas has too little entropy 
for its new position. This means that it is denser and cooler than the gas it replaces, making the radius of this part of 
the star smaller, and reducing the mass transfer rate. To reach thermal equilibrium, this gas must absorb some heat 
from the stellar radiation field, on a thermal timescale. In doing so it expands, returning this part of the star to the 
equilibrium radius for its new (slightly smaller) masq/4| 

The argument above shows that mass transfer from a more massive radiative companion cannot happen on a 
timescale shorter than the star’s thermal timescale, fy, ~ GM; [RL (where R3, Lz are the star’s radius and luminosity 
respectively). But this is usually much shorter than the nuclear timescale in stars of mass a few Mo, and we find 
values -M> ~ M3/tn as large as ~ 1077 — 1074 Mo yr“, depending on the stellar mass M3. 

Rates like this are far higher than the value Mgaa defined above, which would give the Eddington luminosity if it 
could be accreted and so it is very likely that the super-Eddington excess is likely to be expelled. 

As we mentioned in the Introduction and Section 3.6] the neutron star in the long—period (P = 9.84 d) low—mass 
X-ray binary Cyg X-2, appears to have expelled most of the mass transferred to it in this way. The distinctive feature 
of this system is the unusual state of the companion star. Combining spectroscopic information, photometry of the 
ellipsoidal variation of the companion, and the geometry of the Roche potential, showed that 
this star had a low mass M2 = 0.5-0.7 Mo, and yet a large radius (R2 ~ 7 Ro) and very high luminosity (L2 = 150 Lo). 
This is just what one expects at the end of an episode of thermal—timescale mass loss: the hot central region of what 
was initially a much more massive star has been exposed, and has had no time to cool. The evolution followed by Cyg 
X-2 is called early massive Case B, in the terminology of{Kippenhahn and Weigert| (1967). 

The most important feature of this evolution is that the neutron star accretor has evidently gained almost none of 
the ~ 3 Mo transferred to it by the companion: its current mass is only 1.78 + 0.3 Mo. This expulsion cannot have 
come from common-envelope evolution — the orbital binding energy released in shrinking the binary to the separation 
given by the current period of almost 10 days is far too small to expel any matter. This leaves only accretion energy 
— the accretion disc around the neutron star must have expelled most of the transferred mass. This is just the process 
discussed by [Shakura and Sunyaev|(1973), and summarized in Section[3.3] As we discuss in Sections [3.8]and[3.14] it 
leads to the production of a quasi-spherical wind with typical outflow velocity v ~ 0.1c, with open funnels along the 
disc axes which collimate a large fraction of the accretion energy released as radiation. In this picture, ULXs are the 
evolutionary stage following the wind—fed HMXB phase, once the massive companion star fills its Roche lobe. This 
seems particularly likely for the pulsing PULXs, which contain neutron stars very similar to those in HMXBs. 

Although thermal—timescale mass transfer is very probably the cause of ULX behaviour in many cases, it is not the 
only way to make a ULX. Evidently any process that produces a significantly super—Eddington mass transfer rate can 


'4Tf the star instead has a convective envelope, and so constant specific entropy, mass transfer can reach very high (near—dynamical) values, and 
the system develops a deep common envelope. This may lead to a merger of the two a King and Begelman| (1999) 1999) show that this is very 
unlikely to happen with a radiative envelope, even for transfer rates as high as ~ 1074 Mo yr~ 
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do this. We have already mentioned that some Be — X-ray binary systems become ULXs from time to time, because 
of dynamical effects on the Be-star disc. Three other possible ways involve mass transfer on the nuclear timescale 
of a massive donor (Rappaport et al. (2005), transient outbursts involving the thermal—viscous disc instability 
Hameury and Lasota| (2020| see Sect. [3.16), and gravitational—wave driven mass transfer from a white dwarf 
onto a neutron star (King\|2011) — a so—called ultracompact binary (UCB). UCBs are often found in globular clusters 
(see [Dage et al.| and references therein), where they are formed by dynamical capture, so these sources offer 
a possible explanation for ULXs sometimes claimed in elliptical galaxies. We will discuss other evolutionary paths, 
especially for weak (Ly = 3 x 10°’ erg s7!) ULXs in Section[3.13} 


3.12. ‘Super-Eddington Accretion’ 


These possibilities illustrate that ULXs are likely to be binaries containing a compact accretor in an epoch of 
strongly super-Eddington mass transfer. It is worth emphasizing that the commonly used phrase ‘super—Eddington 
accretion’ is highly ambiguous, and a frequent cause of confusion. As the example of Cyg X—2 shows, the accretors 
in ULXs do not gain mass at significantly super—Eddington rates: it is the mass supply rate which is super—Eddington, 
not the accretion. Accordingly, the phrase ‘super—Eddington mass supply’ is strongly preferable. 


3.13. Population syntheses 


Population synthesis of ULXs was first performed by|Rappaport et al.|(2005) who showed that the population of 
ULXs in spiral galaxies can be explained by short phases of high mass-—transfer (with nuclear or thermal timescales) 


in black-hole plus evolved-star binaries. However, they neglected the pre-supernova evolution phase in their cal- 
culations and did not take into account other types of binaries, such as systems containing neutron-star accretors, 
or evolved donors. Their study was expanded by who completed it by predicting the 
observational properties of the ULX population. (2010) used the StarTrack population synthesis code 
to show that the bulk of ULXs can be interpreted as the high-luminosity tail of high- 
mass X-ray binaries. Neutron-star accretors were included in population studies of ULXs only after the discovery of 
PULXs. 

performed a proof-of-concept study where they showed that any ULX, including the 
most luminous ones, and neutron-star systems, can be a short-lived phase in the life of a binary star. The detection 
of double compact object mergers (Abbott et al.||2016) triggered investigations of the connection between binary 
compact objects and ULXs (e.g., [Marchant et al.|/2017}/Klencki et al.|/2018). 

A problem encountered when performing ULX population synthesis is the unfortunate choice of the threshold 
ULX luminosity at 10°° erg s~!, because it fails to separate “standard” X-ray binaries (XRBs) from “generic” ULXs. 
At least five transient Galactic XRBs reach luminosities > 10°? erg s7! (2016). 
suggested that one should take 3 x 10°’ ergs~! as the minimum luminosity defining ULXs whose nature is 
“contentious”. The idea was to ensure that if the accretor was a black hole (BH), the apparent luminosity would be 
super-Eddington. But we now know that stellar evolution can produce black holes with masses > 30 Mo 


2010} {Abbott et al.|/2016), up to 50 Mo under favorable conditions (see, e.g., /Belczynski|2020} the fact that 


the maximum value of black—hole masses observed in X-ray binaries is ~ 20 Mo follows from the limited range of 
metallicities scanned by X-ray observations [Olejak et al.[2020). So only sources with luminosities > 6 x 10° erg s~! 
qualify as (apparently) super-Eddington for a BH resulting from stellar evolution. BH) | The current outdated but 
default definition of ULXs as L > 10% ergs~! tends to coerce population syntheses to use it, at the cost of producing 
“excesses” of sources with luminosities which in most cases are actually sub—Eddington. 

Existing population syntheses do not include Be-X systems, and until recently they considered only Roche-lobe 
overflowing (RLOF) companions of the compact accretor, so excluding wind—fed systems. The second omission was 
corrected by (2021), and wind-accreting ULXs will be discussed below, but there is not much 
hope that the difficulties inherent in the first (Be-X) will be overcome any time soon. Although current models of Be 
stars allow simple modelling of their discs at radii far from the stellar surface, they also contain hidden assumptions 


that are not physically plausible (Nixon and Pringle} (2020). If the model for Be—star discs proposed by the latter 


'The limit can be even higher for He-rich accretion, and for the most massive BHs (e.g.|Belczynski||2020 
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authors (the disc material originating from small-scale magnetic flaring events on the stellar surface) is correct, it will 
be problematic to incorporate it in a population—synthesis code. 


(2017|/2019) performed comprehensive simulations of ULX populations by using the StarTrack 
population synthesis code (Belczynski et al. with significant updates described in/Dominik et al.|(2012); 
(2014). They simulated the evolution of 2 x 10’ binary systems for every model and scaled the 
results to a Milky-Way equivalent galaxy (Mmwec = 6 x 10!° Mo; [Licquia and Newman| 2015). Two star formation 


cases are considered: at a constant rate of 6.0 Mo yr“! for 10Gyr, and 600 Mo yr™! in bursts of star—formation lasting 
100 Myr). They consider three metallicity prescriptions: solar (Zo), 10% of solar ( Zo/10), and 1% of solar ( Zo/100). 
The fiducial accretion model for ULX evolution in this study assumes the SS73 windy solution and beaming described 
by Eq. (55), with a saturation limit at 7 = 150. For higher mass accretion rates, the beaming is assumed to be constant 
and equal to b ~ 3.2 x 10-3 (0 ~ 9°). For comparison, seven other models with different accretion and/or beaming 
configurations are considered. In general, the results do not depend strongly on the assumptions about the accretion 
mode or the beaning model, except for the case where beaming is assumed to be constant (b = 0.1) which, by lowering 
the ULX luminosity threshold, results in neutron—star ULX always dominating the population at late times. 

It is known from observations, that ULXs are often associated with low metallicity environments (e.g. 


2015} Massive black holes apparently form 
more easily in low—metallicity environments (e.g., 2009; Belczynski et al. 
2010), and RLOF onto a compact object is therefore more frequent than in high metallicity cases (Linden et al.||2010). 
(2017) found that metallicity has a strong impact on the number of ULXs, but only in star— 
forming regions. They show that the number of black—hole ULXs increases significantly for Zọ/10 in comparison 
with Zo. But the number of neutron—star ULXs is virtually unaffected by metallicity. Interestingly, for models with 
the lowest investigated metallicity (Zo/100) they find fewer ULXs than for Z = Zo/10. This suggests that the relation 


between the number of ULXs and metallicity is not monotonic (see |Prestwich et al.||2013). The reasons for this 
inversion are currently not understood. 


The general conclusions of the ULX population synthesis by (2017) are: 


e ULX with NS accretors dominate the post-burst ULX populations, and those with constant stellar formation 
(duration > 1 Gyr) in high-Z environments, consistent with current understanding of binary evolution. Neutron— 
star ULXs are present in significant numbers (= 10%) also during the star—-formation bursts and in lower-Z ULX 
populations. 


e ULXs appear in a very specific sequence after the start of the star formation (t denotes the time since the 
beginning of star formation): 
1. t= 4-40 Myr BH-MS (5.6-11 Mo), 
2. t = 6-800 Myr NS -MS (0.9-1.5 Mo), 
3. t x 430-1100 Myr NS-HG (0.6-1.0Mo), 
4. t x 540-4400 Myr NS-RG (~ 1.0 Mo); 


(BH-black hole, NS — neutron star, MS — main sequence, HG — H-rich Hertzsprung gap, RG — red giant.) 
e Neutron—star ULXs may reach luminosities as high as those of black—hole ULXs (Lx max > 10% erg s71). 


e The most luminous ULXs (Lx Z 104! erg s7!) contain HG donors (blackhole ULXs; Mz ~ 1.2-3.7 Mo) or 
evolved helium stars (neutron—star ULXs; M ~ 1.7-2.6 Mo), which overfill their Roche lobe and transfer 
mass on a thermal timescale. They form typically within 15-75 Myr after the ZAMS. 


Finally, as seen in Table [6] the ULX population is in all cases strongly dominated by sources with Ly < 3 x 
10% erg s™!. 

The results of (2017) describe the evolution of the intrinsic ULX population, but since at least 
these sources are supposed to be beamed, the observed ULX population is different. (2019) use 


the results of the population synthesis described above to derive the properties of the observed ULX population. They 


use the fiducial model defined in|Wiktorowicz et al.| (2017) but remove the beaming saturation—condition because 
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Table 6 Number of ULXs per MWEG 100 Myr after the start of stellar formation 


Metallicity > 10% ergs! >3x10%ergs! >10%ergs! > 10% ergs! 
#ULX 407 65 14 0.5 

Zo #BHULX 372 40 4.7 0.38 
#NSULX 35 25 9.4 0.13 
#ULX 7281 1940 200 12 

Zo/10 #BHULX 7200 1900 180 11 
#NSULX 81 40 12 0.43 
#ULX 5000 2400 390 15 

Zo/100 #BHULX 4900 2400 380 15 
#NSULX 16 12 9.0 0.075 


(Adapted from 2017). ‘MWEG’ = ‘Milky Way Equivalent Galaxy’. 


extremely beamed sources are not only hard to observe, but also extremely rare and short-lived so have no influence 
on the final results. 

discuss the full range of possible ways of making a ULX. This paper also estimates that 
at low redshift, about 50 per cent of merging BH—BH progenitor binaries evolved through a ULX phase. They find 
that observed ULXs with black—hole accretors typically emit isotropically (b = 1) and undergo nuclear—timescale 
mass transfer, whereas those with neutron-star accretors are predominantly beamed (typically b = 0.7 — 0.2) and 
in most cases the mass transfer occurs on a thermal timescale. They also show that beaming depends on the stel- 
lar environment; very young (burst) populations (age < 10 Myr), dominated by black-hole ULXs, are significantly 
beamed, while black—hole ULXs in older stellar populations are usually isotropic emitters. In contrast, the majority 
of neutron-star ULXs are always beamed, irrespective of the stellar environment. The ratio of neutron—star ULXs to 
black—hole ULXs is higher in the total sample than in the observed sample. For continuous star formation, black—hole 
ULXs typically outnumber the NS ULXs in the observed sample. Black—hole ULXs also outnumber neutron-star 
ULXs in the observed sample for burst stellar formation at early times, but after the star formation burst, neutron star 
ULXs tend to dominate the observed population. Here the observed neutron—star ULXs represent only 20% of the 
total neutron—star ULX population, and many are expected to be obscured (in the absence of precession, which may 
act to bring some into view, see, [Dauser et al.|[2017} {Middleton et al.|{2018). 

(2019b) pointed out that the growing number of ULXs with (tentatively identified) red—supergiant 
donors (for which RLOF is dynamically unstable) motivates the study of wind-driven mass-transfer, previously ne- 


glected in the context of ULXs. Including wind accreting X-ray binary models into StarTrack, 
(2021) showed that wind-fed ULXs can constitute a significant fraction of all ULXs, and in some environments may 


even be the majority. Bondi-Hoyle-Lyttleton accretion is a standard form of accretion in StarTrack (Belczynski 
and they additionally adapt a wind RLOF (WRLOF) scheme hom e A 
(2007), already used by [fkiewicz et al.] (2019), to analyse the relation between SNIa and wide symbiotic binaries. 
With somewhat optimistic assumptions for wind—accretion efficiency, they suggest that wind-fed ULXs should not be 
neglected. The wind—fed ULX sample they found contains a significant fraction of RSG companions, supporting the 
suggestion that the apparent superposition of some ULXs and RSGs results from coexistence as a binary. Although 
some of these systems will evolve into double compact objects, none of the ULX systems with a RSG donor is a 
viable progenitor of double compact object mergers with fmerge < 1OGyr, because they have large orbital separations. 

It is important to note that models of wind accretion and emission suffer from serious uncertainties. This problem 
cannot be solved only by extensive numerical simulations, but also requires systematic quantitative observational data. 
and [Wiktorowicz et al.| (2021) do not consider the origin or evolution of neutron- 
star magnetic fields, considering that this problem contains too many unknowns to be worth including in population 
synthesis. 


(2020) do not share this view, and perform population synthesis of binaries containing magnetized 
neutron stars only. They assume solar abundances and neglect magnetic—field evolution. They consider two accretion— 


flow models: the standard|Shakura and Sunyaev| (1973) model adapted to magnetized accretors by |Chashkina et al. 
(2017) and for the supercritical accretion they use the|Chashkina et al.|(2019) “advective” model (for the suitability of 
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advection models to describe accretion onto neutron stars see Sect. (3.8.2). The results are presented as a P, —L, figure 
on which the observed values of 10 PULXs are superposed. Although Be/X stars are not present in the modelled 
population, the authors use also four Be/ULX data points. When the data points fall into the regions where the 
calculations predict the existence of PULXs, they consider this as a success of their model. Which is the case for 
all accretion modes, only NGC 5097 ULX-1 shows, as usual, some resistance. They conclude that that standard 
evolution of close—binary systems and accretion onto magnetised neutron—stars can quantitatively explain the observed 
properties of PULX i.e. their X-ray luminosities, spin periods, orbital periods and masses of visual components, 
without assuming beamed X-ray emission. In a model galaxy with star formation rate 3-5 Mo yr! “there can exist” 
several PULXs. If PULXs are not beamed, it is not clear why they are invisible in the Milky Way. (The single PULX 
observed in the Galaxy is in a Be-X binary — because of their uncertain evolutionary status, such systems are not taken 
into account in population syntheses.) assert that discovery of powerful winds from PULXs with 
Lx ~ 104! ergs”! may be a signature of super-Eddington accretion onto magnetised neutron stars. If so, this raises 
two questions. First, in the absence of beaming why would such sources not be observable as bright X-ray sources? 
Second, with such powerful outflows, how one can avoid beaming? Finally, although these authors say that, in their 
models, the PULXs represent “a subset of neutron stars at the stage of disc accretion in close binary systems with 
luminosity Ly > 10% ergs~!”, they take no account of the other distinctive parameter that singles out the PULXs: 
their extremely high spin frequency derivatives v (see Fig. [13). 


3.14. Winds and Feedback from ULXs 


We have seen (in Section [3.6) that the characteristic feature of the disc-wind beaming picture of ULXs is that 
their accretion discs eject the super-Eddington part of the mass transferred to them. If -M >> Mgaa this is of course 
most of the transferred mass. In all cases we expect a quasispherical accretion disc wind, with most of the emitted 
radiation escaping through open funnels along the disc axis. The dynamics of the wind are controlled by radiation 
pressure, and therefore by the electron—scattering optical depth distribution of the outflow. For modest Eddington 
factors ri = M/ Mgaa = 1 we expect that the photons produced by the conversion of accretion energy will find the 
open funnels around the rotational axis of the disc and so escape after only a small number of scatterings. The front- 
back symmetry of non-relativistic electron scattering means that, on average, a photon gives up all of its momentum 
(but not energy) in each scattering event. Then we expect that the accretion disc wind should have an outgoing 
momentum of order that in the original radiation field, i.e. just the Eddington momentum 


f L 
My ~ E (105) 
c 
so that the outflow velocity is 
v= ter Ole, (106) 


where 7 = Leaa/ Megaac? = 0.1 is the radiative efficiency of accretion (King and Pounds} |2003). Winds with the 
properties (105106) are seen in many AGN (see|King and Pounds}\|2015} for a review), where the Eddington factor 


is unlikely to be large. 

But in ULXs, the condition 7 >> 1 is a distinct possibility, and here we expect considerably more scatterings 
before a photon finds the open funnels around the disc axis and escapes. Instead we expect that the radiation field may 
put a significant fraction of its energy into the wind, i.e. 


1. 
zv” = l’ Lead, (107) 


with /’ ~ 1. This now leads to an outflow velocity 


»\1/2 
Ha (7) 7 (108) 


m 


(King and Muldrew{|2016). Despite the different physics involved, the two types of winds have surprisingly similar 


velocities, 1.e. 
v, v = 0.05c — 0.2c (109) 
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for values 10 < m < 100. Observations show the presence of winds with velocities of this order in most (~ 70%) 


ULXs (Sect. |2.2.1). 


An observer not located very close to the central accretion disc axis (i.e. not in the narrow beam defining the ULX) 


would observe the photosphere of the quasi-spherical wind. |King and Muldrew}(2016) show that this has radius 


Ron = 104 x 1? Mono Mio km (110) 
and effective temperature 
Tee = 1 x 1067 4M aM, K (111) 


where M,,19 is the wind mass outflow rate in units of 10r, 70,1 is the accretor efficiency in units of 0.1, and Mjo is the 
accretor mass in units of 10 Mo. 

Objects almost precisely like this — blackbody emission at temperature ~ 100 eV from photospheres of order 
(110) are observed, in the form of the ultraluminous supersoft sources (ULSS, [Kong and Di Stefano] |2003). Evidently 
these are ULXs viewed ‘from the side’. The well-know extreme Galactic source S8433 is probably a member of this 
class, as we discuss below. 

ULX/ULS winds are quasispherical, and hypersonic with respect to the interstellar medium (ISM) surrounding 
the ULX, and so must be slowed in a strong reverse shock against it. A weaker forward shock propagates into the 
surrounding ISM ahead of the contact discontinuity between the ULX wind and the ISM, sweeping this gas outwards 
in a shell. As in the similar problem for supermassive black holes (see[King and Pounds} |2015] for a review) the nature 
of ULX feedback on the interstellar gas depends on whether the shock slowing the ULX wind is cooled by radiation 
losses within a flow timescale or not. In the first case, where the shock is efficiently cooled, most of the energy of the 
wind is lost to radiation, and the wind transmits only its momentum to the ISM. If instead shock cooling is inefficient, 
the outflow gives all of its energy adiabatically to this gas. These two cases are often called momentum-driven and 
energy—driven respectively] 

For ULXs, we will see that feedback is in practice always energy—driven, which in turn means that ULXs generally 
have a significant effect on their surroundings. The temperature of the reverse shock slowing the ULX wind is 
T, ~ myu?/k ~ 10!° — 10!! K, where u ~ v,v’, and k is Boltzmann’s constant. Since this is far higher than the 
T < 10° K radiation field of the ULX, inverse Compton scattering off the ULX radiation field is the main cooling 
mechanism, as in the SMBH case. (The thermal bremsstrahlung timescale is far longer.) King (2003) shows that the 
ratio of Compton cooling time and flow timescale R/v is 


2 
te 2{me\ R 
Pe : (112) 
thow 3 Mp Rg v 
at distance R from the ULX. From (109) this defines a cooling radius 
Re = 5x 107R; ~ 5 x 10 cm, (113) 


which is generally smaller than the ULX binary’s orbit. Shock cooling is never important for ULX outflows, and all 
ULX feedback is energy—driven. 

We can easily see the effect of this feedback. If the interstellar gas has uniform density p and the shock exerts 
pressure P on it, momentum conservation at the radius R of the interface between the ULX wind and the ISM gives 


d [4r : 
4nR?P = — | — RoR 114 
7 s$ ph) (114) 
so that 
P Tos o 
— = —-RR +R’. (115) 
p 3 


16In SMBH accretion, feedback switches from momentum- to energy—driven once the SMBH mass reaches the M — e value M œ o*. This fixes 
the limit to SMBH mass growth, since the ISM gas ultimately fuelling its growth is expelled at this point, in a galaxy—wide high-speed outflow. 
See King & Pounds (2015) 
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Since energy is transmitted adiabatically from the wind we have 


d d (4nR° 
— |2nR°P| = Lea - P— 
ay 2R P| = Lea Al 3 ) 


Eliminating P gives the equation of motion for the shock as 


L E os 
pid WR 4 SRR +5R°R3 (116) 


270 


All solutions of this equation tend to the attractor 


125 Lea ais 
R =| —— t 117 
1012p í an 
(cf Eq. |18) or numerically 
1/5 
L 
R= so(=2) £) pe, (118) 
P-25 


where L39 = Lpga/10*? erg s7}, ts = t/10° yr, and p-25 = p/ 10-25 gcm™. The conditions correspond to a ULX with 
a lifetime of order the thermal timescale of a massive companion star, driving a wind into an ISM of number density 
~ 0.1 cm™ H atoms per cm?. The shell of swept-up ISM gas expands at a speed 


, L39 e -2/5 -1 

k = 300( =] t ~ kms, (119) 
P-25 

so the forward shock into the ISM is much weaker than the reverse shock slowing the ULX wind (cf Eq. |109). For 

longer thermal timescales or more powerful wind mechanical luminosities, L39 ~ 10, and the predicted nebulae can 

have sizes ~ 500 pc. Observations (Sec. show that many ULXs are associated with wind—blown nebulae of order 

50 - 500 pc in size. These are significantly larger than supernova remnants (~ 3 pc). 


3.14.1. SS433 

The well-known extreme Galactic binary system SS433 fits closely into this picture. show 
that it probably has a mass transfer rate implying 7,19 ~ 10°, and ts ~ 1. The mechanical luminosity is then. 
~ 10°’ ergs™!. For a 10 Mo black hole accretor, Eq. implies a spherical nebula (excluding the ‘ears’ formed 
by the precessing jets) comparable with the observed ~ 30 pc size. note that this implies a 


quasispherical wind speed close to the observed velocity width ~ 1500 kms ~ of the so-called stationary H alpha line. 
Despite its very high mass transfer rate, SS433 is a remarkably weak X-ray source —{Watson et al.| (1986) estimate 
Lx ~ 10% ergs™! and a mechanical luminosity ~ 10% erg s™!, and it is likely that the softer (< 10 keV) X-rays come 
only from the prominent jets the system produces (whilst the harder X-ray emission may originate from scattering 
inside the wind-cone, [Middleton et al./202Tb). In a detailed study, show that SS433 has a 
viewing angle appropriate to a ULS, but is not observable as one: from (111) we have kTe ~ 10 ~,eV, which is too 
soft to be detectable as this system is heavily reddened (Ay ~ 7). SS433 is distinctive because of the huge periodic 
redshift and blueshifts of the H—alpha lines emitted by its precessing jets. These evidently result from its extreme 
Eddington factor, and precess (and are baryon—loaded and slowed to 0.26c) probably because they are deflected by a 
thick precessing central disc. As a mark of this, the precession rate of the jets does show power at the binary orbital 


period (Begelman et al.||2006). 


3.15. Be—star ULXs 


The luminosity of ULXs is generally variable, but two classes of ultraluminous X-ray sources are transient. One 
of these consists of systems appearing in the transient X-ray sky only as ultraluminous sources. We discuss these in 
the next section. 

The other class is a small set of systems which normally appear as high—-mass X-ray binaries, but occasionally 
become ULXs. This class has only recently been recognised, but actually includes A0538-66, the first system observed 
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to have an apparently super-Eddington X-ray luminosity (White and Carpenter)| 1978), and also the Milky Way’s only 
known ULX, Swift J0243.6+6124, first detected on 3 October 2017 during a giant X-ray outburst 2017 


see Table [2). In all these systems, a neutron star, often with a noticeable magnetic field, accretes from a much more 
massive Be (or Oe) star in a wide eccentric orbit. 

Be stars are unusual in having a large disc of material surrounding their rotational equator. The origin of the disc 
is unclear, but is evidently related to the star’s rotation, and possibly mediated by small—scale surface magnetic fields 
for a recent discussion). Be X-ray binaries generally produce most of their X-ray emission 
when the neutron star is near pericentre and interacts with the circumstellar disc. The X-rays therefore generally 
indicate the orbital period. 

But it is now known that, on a much longer timescale, this interaction is much stronger, and the system becomes a 
ULX. The origin of this behaviour is probably von Zeipel—Lidov—Kozai (ZLK) oscillations of the Be star disc 
fet al.|/2014). Highly misaligned test particle orbits around one component of a binary cyclically exchange inclination 
for eccentricity |1962). This process conserves the component of the angular 


momentum orthogonal to the binary orbital plane, i.e. 
q- ee cos iy = constant. (120) 


Here i, is the inclination of the particle orbit to the binary orbit and e, is the test—particle eccentricity. A test particle 
initially on a circular strongly misaligned orbit oscillates to closer alignment. Then | cos i,| increases, requiring larger 
ep. The oscillations occur only if the initial inclination i, inclination satisfies 


cos? ipo < z: (121) 
so 39° < ipo < 141°, and the maximum eccentricity that an initially circular test orbit can reach is 
5 1/2 
Emax = ( =a cos? in] : (122) 


(2014) show that this test-particle behaviour extends to full hydrodynamical discs, which undergo 
global KL cycles of the disc inclination and eccentricity. The spins of Be stars are known in general to be strongly 


misaligned from the orbits of their neutron star companions in Be — X-ray binaries, probably as a result of the 
asymmetric supernova explosion producing the neutron star. At maximum disc eccentricity, it is likely that the neutron 
star accretes much more gas from the Be-star disc, and this may turn it into a ULX. Global KL cycles typically have 
periods an order of magnitude longer than the binary orbit, so the transition to the ULX state is relatively infrequent. 


3.16. Modelling transient ULXs 


Accretion discs in low-mass X-ray binaries are subject to a thermal—viscous instability when the mass transfer 
rate from the companion falls below a critical value depending mainly on the size of the disc (orbital period; 


Lasotal 2001 
(van Paradijs 
1996) and the outburst development (King and Ritter||1998} 1999}|2001). 


The thermal—viscous disc instability is naturally associated with low mass transfer rates, but since the critical 
accretion rate increases strongly with the disc size, for large enough binary systems, even super—Eddington mass 


transfer can correspond to an unstable disc 2015| /Hameury and Lasotal |2020). So for sufficiently 


large orbital periods, accreting systems have large unstable discs and can produce outbursts bright enough to explain 


some ULXs (King et al.||2000b). Although as discussed in Section it is very difficult to determine the binary 


parameters of these distant systems, it is clear that many of them have giant or supergiant companion stars (Section 


2.5.1) and have orbital periods longer than 2 days. Several ULXs are classified as transient (see e.g. |Brightman et al. 
2020a). Usually their light curves are only sparsely covered, but the recent well sampled observations of a transient 


ULX source in the galaxy M 51 by|Brightman et al.|(2020a) provides a light curve that can be compared with model 


predictions. Some long-period, sub-Eddington X-ray transients show long-lasting outbursts (tens of years), so some 
apparently steady ULXs could actually be transients observed during long outbursts. 
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The critical rate below which a disc is unstable is well approximated by: 


. R 2.1 
944-04 p-0. out 2 
Mt, ~ 24x 10M; n “| es, (123) 


where Rout is the outer disc radius and fir = 7C/(5 xX 1074) and C is defined through 
mMc? 

4nr? ` 

C contains all the physics of the irradiation process, and m. (2001) found that the light curves of low- 
mass X-ray transients are reasonably well reproduced, and (2012) found that the corresponding stability 
criterion provides the observed division of sources into steady and outbursting, if one uses a constant 1C of order 


10-3. The radiative accretion efficiency 7, is defined as the ratio of the true luminosity (without the beaming factor) to 
c?, and is therefore: 


ot =C (124) 


m = 0.1 if m< l, 
=0.1(1 +lnň)/  ifm>1. (125) 


From Eq. it follows that neutron-star ULXs may be unstable for super—Eddington mass transfer rates at orbital 
periods = 10 days (slightly shorter for black-hole ULXs). But if the front fails to reach the outer disc edge, the 
peak outburst luminosity is given by Eq. with Rout replaced by Remax < Rout, where Ry,max is the maximum 
radius reached by the instability-triggering heating front. In this case, the peak luminosity can be up to 50 times 


larger (Hameury and Lasota}|2020). This means that wide binaries in which the mass transfer is sub-Eddington might 


become transient ULXs at, and close to, the maximum outburst luminosity if their discs are unstable. 


Hameury and Lasota| (2020) derived analytical formulae describing typical X-ray transient light curves with great 


accuracy. They apply also to transient ULXs. When the heating front brings the whole disc to a hot state (Ri max * 


Rou), the M(t) relation during the decay from outburst maximum is very well described by (Ritter and King}/2001 17 


1713 
M = Mafi, (126) 
0 
where the integration constant tọ is given by: 

to = 3.19099 MARY MIN, yr, (127) 

with Mmnax.i9 = Mmax/10!° gs7!. Here to is not the characteristic timescale of disc evolution, which is instead 
_ Ma _ O.81y F-93037 70-15 0.62 9-08 128 
C= ON 1 Sie Riz? %2 Y5 (128) 


where M4 is the disc mass and f = M max /M*,(rou) > 1 (for the derivation see|/Hameury and Lasota}|2020). One can 


then relate the maximum outburst accretion rate to the characteristic decay time. 


irr 


3.39 
Minax = 2.0 X 101922! £2.02 m7165 (5) -1.01 z si (129) 


Then when both Mmax and T are known from observations, this relation determines f, and hence Mta and the size of 
the accretion disc. 

This phase of the outburst corresponds to a viscous decay, during which X-ray irradiation prevents the formation 
and propagation of a cooling front. This ends when the decreasing accretion rate falls to the critical value M+. and 


crit 
the front starts propagating into the disc. 


17The corresponding formula in the original publication contains a misprint. 
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If the heating front is unable to reach the outer regions — which happens quite often in huge discs — the viscous 
decay phase is missing and, from the start, the decay-from-maximum M(t) relation corresponds to the cooling—front 
propagation, and has the form 


3.39 


B at (130) 


M = 1.31 x 10% a2] M!S fy! e is 


where ¢, is a constant that is determined by the condition that, when the cooling front starts, M is equal to Mt. at the 
maximum transition radius. Then t can be written as: 

r -1 470.37 £0.15 -0.8,.0.62 

= 4.64 € Mi fia 2 ri yi, (131) 
and £ is a constant to be calibrated by numerical simulations; it is found to be in the range 2.47 < € < 8.0. Here 4 


is equal to the duration of this phase of the outburst, and can be directly compared to observations. As before, there 
exists a relation between Mmax and f): 


3.39 
Mmax = 3-1 x 10033! MTS a fao gs”, (132) 
which is essentially the same as in Eq. (129), with f = 1. 
It is important to realise that, for super-Eddington accretion rates, M(t) does not directly represent the observed 
Lx(t) light curve, as the relation between the luminosity and accretion rate is no longer linear. Then when considering 


an outburst whose apparent luminosity is larger than Lggq|Hameury and Lasota| (2020) use the following relations 


mn 
L; = (1 + Inm) bic Leda if m> 1 


= M Leda if m< 1, (133) 


where b = (1 + m/b! is a beaming term; the larger the beaming parameter b, the larger b and hence the smaller 
the beaming effect. Eq. corresponds to b = b/m? with b = 73. Since the beaming factor b/m? applies only for 
in > Vb, in time-dependent calculations, one replaces it by (1 + m/b)! in order to get a smooth transition with the 
case where beaming is negligible. 

Further, when rin > 1 so that the apparent luminosity scales approximately as rn, the decay time of the luminosity 
TL is one half of the decay time of the accretion rate, and Eq. has to be written as: 


3.39 
imax = 132 a27! £2.02 M7265 | = ) Go (134) 
yr 


If not all of the disc reaches the hot state, the total outburst duration is related to the maximum mass accretion rate 
via Eq. (132), and we must have: 


7 \3.39 


Himax = 24.0054 MT [ Sty) fin 


where to now refers to the outburst duration. 


3.16.1. M51 XT-1 

M51 XT-1 is currently the only ULX transient source whose light curve is suitable for comparison with the 
prediction of the disc—instability model. Figure [19] shows how observational data of M51 XT-1 (Brightman et al. 
compare with the predictions of the disc—instability model when the light curve comes from 
combined with Eq. (133) with the beaming parameter b = 73. Two cases are shown. A 1.4 Mo accreting neutron 


star with a disc size 4.8 x 10!! cm; and a 10 Mo accreting black hole, with a disc extending to 4.9 x 10!! cm. As 
can be seen, the agreement is reasonably good, and as acceptable as the original fit with a power law with index —5/3 
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Figure 19 Observed flux from M51 XT-1 as compared with model predictions for a 1.4 Mo neutron star (red curve) 
and a 10 Mo black hole (green curve). Black points correspond to SWIFT data and blue points to Chandra or XMM- 


Newton observations (courtesy M. Brightman, see also|/Brightman et al.||2020a). 


(Brightman et al.||2020a). These data can also be fitted with a 10 Mo accreting black hole, with a disc extending to 
4.9 x 10" cm. 


In the neutron star case, the maximum accretion rate is 6 x 10!9 gs}, while in the black hole case it is 1.5 x 
10% gs-!. Hence the emission from a black-hole system would not be beamed, but for a neutron-star accretor the 
beaming factor would be b = 0.06, as expected for m = 33. The mass transfer rate can be very roughly estimated to 
be of order 1 — 2% of the maximum accretion rate, i.e. 1 — 3 x 10!8 g sl, 


Hameury and Lasota|(2020) showed that taking b = 200 instead of b = 73 in Eq. (133) also provides an acceptable 


fit, but b = 20 is excluded by the Chandra and XMM-Newton points, thus suggesting only moderate beaming. 


3.16.2. HLX-1 in ESO 243-49 

HLX-1 in the galaxy ESO 243-49 attracted a lot of attention as the best IMBH candidate 
among ULXs. It is the brightest ULX candidate known; its luminosity is variable and can exceed 10” erg s~! at 
maximum. Its association with the galaxy ESO 243-49 at a distance of 95 Mpc is fairly well established (Wiersema 
2010} so there is little doubt about its luminosity (see, however, [Lasota et al.|/2015). 
(2020) conclude that a self-consistent scenario is that HLX-1 is a central massive black hole on the low-mass end of 
the mass distribution of nucleated 10? — 10!° Mo galaxies. Its small host galaxy was accreted and tidally stripped by 
ESO 243-49, leaving only the bare nuclear star cluster near the black hole (e. g- Mapelli et al.[2013). Star formation, 
perhaps encouraged by the merger, could have increased the capture rate for the hole to acquire a companion star on 
which it is currently feeding. Then by the definition adopted in this review HLX-1 is not an IMBH candidate, but was 
the central MBH of a small galaxy, which (unusually) has been accreted and disrupted by a larger one. 

In line with this picture, HLX-1’s historical light curve (see e.g. initially showed a promising 
series of almost annual outbursts, only to lose this quasi—periodicity later on. Subsequent to its discovery, a number 
of so—called QPE (‘quasiperiodic eruption’) sources have recently been found from the centres of small galaxies (e.g. 
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2019} |2020 2021). These can be modelled as tidal disruption near—misses (King|2020): 
instead of being fully disrupted by filling its tidal radius near the MBH (a tidal disruption event or TDE) an infalling 


star fills its tidal lobe only at the pericentre of a very eccentric orbit about the MBH. In the small galaxy GSN 069 
the star is probably a low-mass white dwarf (King/2020), and the predicted CNO anomaly was later found 
et al.|2021). It appears that this model works well for all the currently recognised QPE systems 2022). These 
systems all have very short lifetimes as the orbits shrink under gravitational radiation, which is what drives the mass 
transfer. 

Modelling HLX-1 instead as a transient IMBH source undergoing disc outbursts shows that the outburst timescales 
are incompatible with a central mass of > 10* Mo, as would be required to explain the peak luminosity of 10*! erg s7! 


(Hameury and Lasota|2020)and references therein). 


3.17. Beaming or not? 


The arguments in favour of luminosity beaming in PULXs presented above are simple and fundamental. The 
main one is that the observed X-ray luminosities of magnetic ULXs lie between 6 and 500 Lgaa, with the classic 
(non BeX) systems having Ly > 20Lgaa. Even if mass is supplied at a super—Eddington rate, the luminosity is only 
L = Lgag(1 + Ini), so even 1 ~ 1000 (> 1075 Moyr™!) produces only ~ 8 Lead. 

Explaining these huge excess luminosities by appealing to super—strong magnetic fields appears unpromising. We 
have argued above that the presence of magnetars in PULXs is extremely unlikely. Among other problems, it requires 
a cosmic conspiracy confining magnetars to just these binary systems. Be-PULXs spend most of the time as standard 
pulsed X-ray sources, and clearly do not have magnetar-strength fields. Assertions that the “apparent luminosity [of 


PULXs] is close to the actual one” (Mushtukov et al.|/2020) are rather surprising, especially when it is not explained 


how this is supposed to be achieved. 

refer to a table in{King and Lasota|{2020) to claim that these authors require a = 1/b ~ 20 
for M82 ULX-2, NGC 7793 P13 and NGC 300 ULX1 and ~ 100 for NGC 5. But Table [3] shows that the values of 
1/b for the first three sources are respectively 17, 6, 6, to which one can add 10 for ULX-7 in M51. 


3.18. Measuring mass transfer via the orbital period derivative? 


The most important parameter in understanding ULXs, as in any accreting system, is the instantaneous mass 
transfer rate —M> in the host binary system. Unfortunately, precisely because we do not yet have a clear picture 
of how exactly the mass transfer rate specifies the luminosity of ULXs, it is difficult to measure —M) directly from 
observation. But if we assume conservative mass transfer (i.e. conserving both total mass and angular momentum) we 
see from Section that all quantities such as the companion mass M3, the binary separation a, and by Kepler’s 
law, the binary period P œ a*/*, evolve on the mass transfer timescale M>/(—M>) — see Eq. (104). This makes it 
tempting to use the much more easily measured orbital period derivative P to estimate —M>. The recent papers by 
try to use this method for the ULX M82 X-2. 

As we noted above, the problem of measuring -M is common to all interacting binary systems, not simply ULXs. 
There is a long history of efforts to find —M> by using P and assuming conservative evolution, going back at least 
to early studies of cataclysmic variables (where a low—mass main—sequence star transfers mass to a white dwarf in a 
very circular orbit). But all these attempts give uniformly disappointing results — not least because even the sign of P 
often flatly contradicts that expected from conservative evolution. The underlying reason for these difficulties is that 
relations like hold only if the quantities P, —M are averages over a timescale long compared with the time ty 
for the Roche lobe to move one scaleheight within the donor star. This requirement arises because stars have smooth 
density distributions rather than sharp edges. On timescales shorter than fy the instantaneous value of P is always in 
practice dominated by large but short-term effects, with varying signs often in conflict with the long-term rate. This 
point has been made repeatedly in the literature, going back at least to|Pringle|(1975); see also e.g. and 
Frank et al. (Section 4.4). Pringle (1975) concluded that “... an accurate estimate of the rate of mass transfer 
cannot be deduced from a change of binary period”. 


In response to (2021a),|King and Lasota|(2021) estimated ty ~ 1000 yr for M82 X-2. Then obser- 
vations of P over humanly accessible timescales cannot give information about —Mp. |King and Lasota|(2021) further 


pointed out that direct application of the method of (2021a) to three well-studied X-ray transients 
with known values of P (Nova Muscae 1991, XTE J1118+480, A0620-00, see|Gonzdlez Hernandez et al./2017/and 
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references therein) was possible, but not encouraging. These Roche-lobe filling systems are all currently in the faint 
quiescent state. Under the assumption of conservative evolution, their values of P would instead predict that all three 
of them should currently have high mass transfer rates. Their accretion discs should therefore be steady rather than 
transient, putting them stably in bright states — indeed Nova Muscae 1991 would actually be a ULX. The subsequent 


publication (Bachetti et al.|/2021b) does not answer these points. 


4. Conclusions 


Table [7|summarises how various models of ULX behaviour compare with observations. The case for disc—wind 
beaming appears fairly strong, and is so far compatible with a large body of data. Given that the photon—bubble 
model does not appear promising, this is the only model which can straightforwardly apply to both black—hole and 
neutron-star accretors. 

Other models come into consideration if we allow for the idea that ULXs might not be a single homogenous 
population, but instead consist of various different types of systems (as in the example of long and short gamma-ray 
bursts). IMBH clearly cannot account for PULXs or Be—star ULXs, and have not been widely tested in other aspects. 
Unless their masses are > 3 x 10*Mo they cannot account for the striking universal Lsoft « T~* behaviour of ULX 
soft X-ray components with luminosities > 3 x 10°? erg s~!. 


Table 7 Models of ULXs 


model blackbody L œ T~* [winds & nebulae] evolution luminosity fnJL > 3 x 10” ergs" "/pulsing|BeX|BH & NS 
disc-wind beaming yes yes yes yes yes yes | yes yes 
magnetic beaming no ? yes us yes yes | no no 
IMBH needs M > 3x 10f Mo ? ? ? yes no | no no 
magnetars no no fieldstrength? no no yes | no no 
photon bubble no no ? ? no no | no yes 
oriented jets ? ? ? ? no ? no ? 


Observations of Be — X-ray binaries pose a particular problem for models invoking magnetic beaming, as they 
make regular transitions between normal Be X-ray binary states to PULXs, and back again, when their mass supply 
rates are increased or decreased by Kozai—Lidov cycles. None of the Be—star ULXs is observationally distinct from 
other Be X-ray systems when it is in its usual Be X-ray binary state, and in the ULX state these systems have all 
the same properties as other PULXs, including a markedly more rapid spin up rate (cf Fig[13} None of the Be X-ray 
binaries appears to have a very strong magnetic field. This makes it inescapable that very strong dipole magnetic 
fields are not required for making a PULX. Highly—magnetic PULXs might then at most constitute some fraction of 
the PULX population, but there is so far no clinching demonstration that any PULX contains a neutron star with this 
kind of field. Even then there remains the difficulty that this model requires us to believe that strong—field neutron 
stars in binaries are unobservable unless the companion star supplies them with mass at a super—Eddington rate. 
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Appendix A. “Standard” X-ray pulsars 


To illustrate the fact that PULXs differ from standard X-ray pulsars not only in luminosity > 10°? erg s~!, but also 
in their very large spin-up rates v > 107!9s~?, we used the XRP sample with known luminosities and spin-up rates 


(Table|Appendix A) as compiled by (1985). 
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Table A.8 Observed properties of ‘classical’ X-ray pulsars 


Name Lx(max) [erg s71] P, [s] v [s-?] 

SMC X-1 50x105 [0.71 12.6 x 107" 
Cen X-3 5.0x 10°77 14.84 11.9 x 10-2 
GX 1+4 4.0x 10%7 122 15.2 x 10°? 
4U 0115+63} 3.4x 10°” = 13.61 12.2 x107? 
A 0535+26 1.0x 107 J104 13.0 x107? 
Her X-1 1.0x 107 [1.24 18.5 x 107" 
Vela X-1 1.5x 10% 1283 [5.6 x 107" 
GX 301-2 1.0x 10% [696 13.8 x 10-8 
X-Per 3.9x 103 1835 16.4 x 1071 
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